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Abstract 

Ionospheric  scintillation  of  GPS  signals  threatens  navigation  and  military  op¬ 
erations  by  degrading  performance  or  making  GPS  unavailable.  Scintillation  is  par¬ 
ticularly  active,  although  not  limited  to,  a  belt  encircling  the  earth  within  ±20°  of 
the  geomagnetic  equator.  This  belt  also  hosted  roughly  half  of  the  completed  U.S. 
military  operations  in  the  last  decade.  We  examined  scintillation  data  from  Ascension 
Island,  U.K.  and  Ancon,  Peru,  in  the  Atlantic  &  Americas  longitudinal  sector  at  as 
well  as  data  from  Parepare,  Indonesia  and  Marak  Parak,  Malaysia  in  the  Pacific  lon¬ 
gitudinal  sector.  From  these  data,  we  calculate  percent  probability  of  occurrence  of 
scintillation  at  various  intensities  described  by  the  S4  index.  Additionally,  we  deter¬ 
mine  Dilution  of  Precision  at  one  minute  resolution.  Diurnal,  seasonal  and  solar  cycle 
characteristics  are  examined.  Latitudinal  and  longitudinal  comparisons  are  made. 
Onr  findings  are  consistent  with  previous  research.  Unlike  previous  research,  how¬ 
ever,  we  attempt  to  replace,  or  impute,  missing  S4  values  in  order  to  better  capture 
the  extent  of  scintillation.  In  doing  so,  we  study  data  gaps,  or  holes,  and  characterize 
them. 
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HOLES: 


IONOSPHERIC  SCINTILLATION 
GPS 

AND  IMPUTATION 

I.  Introduction 

“Wherever  you  go,  there  you  are.” 

1 . 1  Overview 

Over  the  centuries  various  tools  and  techniques  have  been  developed  specifically 
to  answer  the  question,  where  are  we?  Accurate  position  information  is  critical  for  a 
myriad  of  human  endeavors,  from  travel  and  trade  to  disaster  relief  and  defense.  This 
fact  has  been  reflected  in  the  refinement  of  navigational  tools  and  techniques  over  the 
centuries.  Increasing  navigational  precision  and  accuracy  has  been  paralleled  by,  and 
to  some  degree  driven  by,  a  concurrent  increase  in  the  lethality  of  weaponry. 

The  United  States  military  has  publicly  promoted  the  precision  of  its  weaponry. 
Phrases  like  “surgical  strike”  and  “smart  bomb”  have  become  part  of  the  popular  lex¬ 
icon.  Understandably,  political  leaders  and  the  general  public  have  increased  expec¬ 
tations,  amplifying  the  political  consequences  of  collateral  casualties  from  a  misplaced 
missile  strike  or  blundered  bombing.  Natural  phenomena  that  degrade  targeting  ac¬ 
curacy,  therefore,  merit  understanding.  This  is  the  thrust  of  our  research.  In  particu¬ 
lar,  we  examined  the  impact  of  a  phenomenon  called  ionospheric  scintillation  on  the 
availability  and  precision  of  information  from  the  Global  Positioning  System  (GPS). 
We  found  that  depending  on  the  temporal  and  spatial  circumstances,  the  impact  of 
scintillation  can  range  from  the  benign  to  catastrophic,  from  virtually  no  discernible 
impact  to  a  total  system  outage  for  the  affected  region. 
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Since  its  inception  in  the  early  1970s,  the  GPS  has  become  the  cornerstone  of 
our  military’s  ability  to  effectively  navigate  and  target.  Twenty  four  GPS  satellites 
orbit  the  Earth  at  at  20,200  km.  Using  signals  from  as  few  as,  but  no  less  than 
four  of  these  satellites,  a  receiver  on  the  ground  can  calculate  the  transit  time  of  the 
electromagnetic  (EM)  signals  from  the  satellites,  and  with  some  adjustment  for  clock- 
error,  generate  position  information.  Today,  GPS  receivers  can  be  found  not  only  in 
weapon  systems,  but  also  in  commercial  vessels,  automobiles,  and  mobile  telephones. 


Before  reaching  a  receiver,  however,  the  GPS  transmissions  must  pass  through 
the  Earth’s  ionosphere.  This  region  of  plasma  envelops  the  Earth  from  beyond  1000 
km  down  to  approximately  60  km  and  is  formed  primarily  through  photo-ionization, 


but  subject  to  other  chemical  and  kinematic  processes  Schunk  and  Nagy ,  2000  .  The 
plasma  refracts  the  GPS  transmission,  delaying  the  arrival  of  the  carrier  wave.  Iono¬ 
spheric  path  delay  occurs  in  a  manner  that  is  reasonably  predicted  and  ameliorated; 
in  fact,  it  can  be  exploited  to  estimate  the  integrated  number  density  of  electrons  in  a 
column  with  a  one  m2  cross  sectional  area  extending  from  the  receiver  to  the  satellite; 
the  Total  Electron  Content  (TEC).  Scintillation  of  the  signal  because  of  ionospheric 
irregularities,  on  the  other  hand,  is  less  predictable  and  more  difficult  to  contend  with. 

Although  scintillation  can  occur  anywhere  over  the  Earth,  certain  regions  are 
conducive  to  the  phenomenon.  One  such  region  encircles  the  Earth  within  ±  20° 
of  the  geomagnetic  equator,  forming  a  “scintillation  belt”  that  encompasses  approx¬ 
imately  1/3  of  the  Earth’s  surface  [ Schunk  and  Nagy,  2000  .  This  region  has  also 
hosted  roughly  half  of  the  completed  United  States  military  operations  in  the  past 
decade,  according  to  data  culled  from  globalsecurity.org.  Not  surprisingly,  receivers 
were  deployed  throughout  the  “scintillation  belt”  beginning  in  the  1990s  to  monitor 


ionospheric  scintillation  at  GPS  frequencies  Groves  et  al. ,  1997  .  The  confluence  of 
these  three  factors  provided  both  the  opportunity  and  the  impetus  to  study  equatorial 
ionospheric  scintillation  and  its  impact  on  GPS. 
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1.2  The  Research 


The  goal  of  our  research  was  to  determine  the  number  and  percentage  of  GPS 
satellites  affected  by  various  levels  of  scintillation.  The  amount  of  error  associated  with 
GPS  measurements  is,  to  a  certain  extent,  a  consequence  of  the  number  and  position 
of  the  available  GPS  satellites,  so  we  also  examined  the  impact  of  scintillation  on  the 
horizontal  dilution  of  precision,  or  HDOP.  We  examined  data  collected  during  solar 
minimum  and  maximum  at  four  stations.  This  work  is  essentially  a  continuation  of 


that  begun  by  Thomas  et  al.  2001  .  That  research  involved  similar  calculations,  with 
the  exception  of  HDOP,  for  stations  in  the  Pacific  during  1998-1999,  prior  to  the  solar 
maximum  in  2001. 


During  severe  scintillation,  receivers  can  lose  the  signal  from  one  or  more  satel¬ 
lites;  civilian  and  military  GPS  users  need  to  know  about  this  phenomenon  because 
navigation  is  threatened.  Statistics  that  do  not  take  into  account  this  phenomena 
risk  underestimating  the  magnitude  of  scintillation.  Assuming  that  data  are  missing 
precisely  because  of  severe  scintillation,  we  could  replace,  or  impute,  the  missing  data 
with  a  proxy  to  indicate  the  severity.  This  is  at  best  a  gamble,  and  at  worst,  a  serious 
distortion.  However,  if  done  carefully  and  the  results  used  with  caution,  it  can  provide 
a  more  complete  picture  of  the  synoptic  behavior  of  this  phenomenon. 

The  imputation  scheme  we  developed  was  able  to  produce  realistic  patterns 
of  scintillation.  Examination  of  imputation  augmented  statistics  from  Ancon,  Peru; 
Ascension  Island,  U.K.;  Parepare,  Indonesia,  and  Marak  Parak,  Malaysia  reflected 
seasonal,  diurnal,  solar-cycle,  and  spatial  dependencies. 

In  order  to  create  a  reasonable  imputation  scheme,  we  had  to  learn  as  much  as 
possible  about  the  missing  data  events,  or  as  we  called  them,  holes  or  outages.  The 
amount  of  missing  data  is  typically  small  and  therefore  ignored.  To  our  knowledge, 
this  is  the  first  research  in  which  time  was  spent  examining  the  outages.  Amplitude 
scintillation  is  typically  characterized  by  the  ratio  of  the  standard  deviation  of  signal 
intensity  to  its  mean,  which  is  called  the  S4  index.  We  found  a  strong  correlation 
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between  the  maximum  value  of  the  S4  index  and  the  variance  of  64  prior  to  an 
outage.  Outages  usually  lasted  less  than  10  min  and  longer  duration  outages  were 
generally  associated  with  low  elevation  satellites.  These  findings  guided  our  design  of 
an  imputation  scheme. 

1.3  This  Work 

The  remainder  of  this  work  is  concerned  with  understanding  the  synoptic  be¬ 
havior  of  scintillation  and  quantifying  its  impacts  on  GPS.  In  Chapter  II,  the  reader 
is  acquainted  with  GPS,  introduced  to  the  ionosphere  and  finally,  provided  with  back¬ 
ground  on  scintillation  and  the  research  to  date.  Chapter  III  is  concerned  with  the 
getting  the  numbers.  We  describe  the  techniques  used  to  discriminate  between  the 
good  and  bad  data,  and  fill  in  the  missing  data.  Here  we  digress  and  examine  the  char¬ 
acteristics  of  holes.  Chapter  IV  contains  the  results.  The  striking  examples  are  here, 
demonstrating  the  temporal  and  spatial  variations  as  well  as  pointing  out  the  utility 
of  imputation.  Finally,  Chapter  V  concludes  with  my  ideas  for  continued  exploration 
of  scintillation  using  ground  based  observations  collected  by  GPS  receivers. 
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II.  Background 


This  chapter  provides  background  on  GPS,  the  ionosphere  and  the  interactions 
leading  to  scintillation.  A  discussion  of  scintillation  and  its  impact  on  electro¬ 
magnetic  waves  follows.  We  end  with  a  survey  of  research  to  date  on  scintillation  and 
its  impacts  on  GPS. 

2.1  The  Global  Positioning  System 

The  first  step  toward  exploring  the  effect  of  scintillation  on  GPS  requires  an 
understanding  of  GPS  itself.  As  mentioned  in  the  preceding  chapter,  a  constellation 
of  24  operational  GPS  satellites  surrounds  the  Earth  at  an  altitude  of  20,200  km,  in 
six  orbital  planes,  with  an  orbital  period  of  12  sidereal  hours.  As  we  imagine  the 
satellites  orbiting  overhead,  several  questions  come  to  mind.  For  instance,  how  are 
satellite  signals  used  to  determine  location?  What  coordinate  systems  are  used?  The 
signals  emitted  by  the  GPS  satellites  are  complex,  carrying  a  great  deal  of  information. 
Which  ones  were  collected  for  this  research?  Additionally,  errors  from  a  variety  of 
sources  accumulate  as  signals  travel  from  the  satellites  to  a  receiver  and  are  processed; 
these  must  be  accounted  for.  The  errors,  in  turn,  influence  the  precision  of  the  location 
estimate.  How  are  these  related?  Finally,  satellite  orbit  data  are  available  in  varying 
degrees  of  accuracy.  Which  will  be  used  for  this  research?  Unless  otherwise  noted, 
the  organization  and  content  of  this  section  was  culled  from  course  notes  provided  by 
J.  Racquet  (EENG  533,  Navigation  Using  the  GPS,  course  handouts,  2006). 

2.1.1  Determining  Location.  One  method  of  determining  location  involves 
transmitting  a  signal  and  recording  the  precise  time  the  signal  was  sent  from  the 
transmitter  and  the  precise  time  it  arrives  at  a  receiver.  With  knowledge  of  the 
signal’s  velocity  through  the  transmission  medium,  the  range  can  be  extracted.  For 
electromagnetic  waves  in  a  vacuum,  ignoring  relativistic  effects  and  assuming  a  fixed 
transmitter  and  receiver,  the  distance  to  the  transmitter  to  the  receiver  is  given  by 
d  —  ct  where  c  is  the  speed  of  light  and  t  is  the  travel  time  for  the  signal  to  reach 
the  receiver.  At  this  point  only  the  distance  from  a  particular  transmitter  is  known, 
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but  not  yet  the  location.  Two  additional  transmitters  are  needed  to  produce  a  unique 
solution  (location).  With  signals  from  all  three  transmitters  providing  the  range 
to  each,  the  relative  position  can  be  determined.  The  intersection  of  three  circles, 
each  centered  on  a  transmitter  with  a  radius  equal  to  the  calculated  range,  provides 
the  position  to  the  transmitters  (for  a  two  dimensional  case).  This  assumes  perfect 
synchronization  of  all  three  transmitter  clocks  and  the  receiver  clock.  This  case  can 
easily  be  extended  to  three  dimensions  by  replacing  the  circles  with  spheres. 


In  practice,  GPS  receiver  clocks  are  less  accurate  and  not  synchronized  with 
the  satellites,  so  clock  error,  5t,  is  introduced.  This  is  the  offset  between  the  receiver 
clock  and  the  satellite  clock  \Hofmann-Wellenhof  et  aL.  2001  .  The  distance  equation 
becomes 

d  =  d  +  cSt  (2-1) 


The  distance  is  typically  referred  to  as  range,  and  range  with  clock  error  is 
called  pseudo-range ,  denote  by  p.  The  consequence  of  clock  error  is  to  require  a 
fourth  satellite  to  accurately  determine  position.  Returning  to  the  three  circles  and 
introducing  a  fourth,  the  position  is  now  at  the  center  of  a  fifth  circle ,  tangent  to  the 
other  four.  The  radius  of  the  fifth  circle  represents  the  Ad  introduced  above.  Again, 
this  can  be  extended  to  three  dimensions. 


Approaching  the  problem  from  another  perspective,  in  order  to  know  location  in 
three  dimensions,  three  coordinates,  x,  y ,  and  z  are  required.  With  three  unknowns, 


three  versions  of  equation  (2.1)  are  needed  to  solve  for  d,  hence  signals  from  three 
satellites.  Having  introduced  a  fourth  unknown  St,  four  equations  (and  consequently 
four  satellites)  are  necessary  in  order  to  produce  a  navigation  solution.  The  resulting 
system  of  equations  can  be  solved  using  methods  described  by  Hofmann-  Wellenhoj 


et  al.  2001,  chap.  9]. 


Until  now  the  discussion  of  pseudo-ranges  to  GPS  satellites  has  taken  place  in  a 
Cartesian  reference  frame  centered  on  the  user.  In  the  long  run  however,  the  typical 
GPS  user  is  not  concerned  with  his  or  her  position  relative  to  a  constellation  of  GPS 
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satellites,  but  rather  his  or  her  position  with  respect  to  Earth  or  locations  on  the 
Earth.  Before  moving  further,  therefore,  three  coordinate  systems  are  introduced  and 
discussed. 


2.1.2  Coordinate  Systems  and  Transformations.  We  are  interested  in  three 
coordinate  systems: 


•  Earth-Centered  Earth-Fixed  (ECEF) 

•  Geodetic  or  Latitude,  Longitude,  Altitude  (LLA) 

•  Local  or  East,  North,  Up  (ENU) 


ECEF  is  a  Cartesian  reference  frame  centered  on,  and  rotating  with,  the  Earth. 
The  LLA  frame  is  based  on  the  World  Geodetic  System  1984  WGS-84  ellipsoid,  a 
global  three-dimensional  coordinate  system  used  with  geospatial  data.  In  the  sim¬ 
plest  terms,  the  Earth  is  represented  by  a  sphere  that’s  ’’squashed”  at  the  poles  and 
bulges  at  the  equator.  The  latitude  (</>),  longitude  (A)  and  altitude  (m),  with  re¬ 
spect  to  the  ellipsoid,  determines  position.  Finally,  the  ENU  frame  is  fixed  to  a  point 
and  oriented  east,  north  and  up.  Later  it  will  be  necessary  to  rotate  coordinates 
from  an  ECEF  to  a  ENU  framework.  The  ECEF  coordinates  must  be  adjusted  with 
respect  to  the  origin  of  the  local  level  frame.  Placing  the  satellites,  as  well  as  the 


receiver  in  a  particular  reference  frame  is  essential  because,  as  noted  by  Hofmann- 


Wellenhof  et  al.  2001,  chap.  4],  ’’For  single  receiver  positioning,  an  orbital  error  is 


highly  correlated  with  positional  error.”  Information  about  a  satellite’s  orbit,  also 
known  as  ephemerides,  is  typically  obtained  from  one  of  three  sources,  almanac  hies, 
broadcast  ephemerides  and  precise  ephemerides.  Ephemerides  allow  a  satellite  to 
be  placed  in  the  desired  reference  frame.  Of  the  three  sources,  almanac  data  are 
the  most  coarse  with  an  uncertainty  of  some  kilometers.  It  is  updated  weekly  at  a 
minimum  and  included  in  the  broadcast  satellite  message  \Hofmann-Wellenhof  et  al 


2001  .  Software  used  for  this  research  ingested  almanac  data  in  YUMA  format  (one 


of  two  almanac  formats)  available  from  the  United  States  Coast  Guard  Navigation 
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Center  at  http://www.navcen.uscg.gov/gps/almanacs.htm.  Additional  information 
regarding  orbit  data  can  be  found  in  the  resources  mentioned  earlier  and  also  the 


encyclopedic  volume  by  Parkinson  et  al.  1996 


On  the  topic  of  orbits,  it  is  useful  to  note  that  the  GPS  ground  track  repeats 
every  sidereal  day,  the  time  it  takes  the  Earth  to  rotate  once  on  its  axis  relative 
to  inertial  space  (the  stars).  The  sidereal  day  is  shorter  than  the  solar  day  by  four 
minutes,  so  a  given  satellite  will  become  available  four  minutes  earlier  each  day.  The 


effects  of  this  behavior  will  become  apparent  on  some  of  the  figures  in  Chapter  IV 


The  alert  reader  may  have  noticed  the  surreptitious  introduction  of  another 
error  source,  orbit  error,  in  the  preceding  discussion.  Rather  than  attempt  to  slip 
any  more  error  sources  in,  we  present  a  brief  summary  with  attention  focused  on  the 
sources  most  relevant  to  our  research. 


2.1.3  Sources  of  Error.  Errors  in  addition  to  those  introduced  by  the  re¬ 
ceiver  and  satellite  can  be  appended  to  the  pseudo-range  equation  [2Tj  which  becomes: 


p  =  r  +  c  (Stu  -  Stsv  +  ^ 


other 


(2.2) 


In  equation  2.2  the  5t  term  from  2.1  was  expanded,  accounting  separately  for  the  re¬ 


ceiver  (user)  error  5tu,  the  satellite  (space  vehicle)  error  Stsv  and  ’’all  others”  Mother- 
The  last  term  is  composed  of  errors  introduced  by: 


•  Multipath 

•  the  Ionosphere 

•  the  Troposphere 

•  Receiver  Noise  and  Resolution  Error 

•  Hardware 


Selective  Availability 


The  italicized  items  indicate  the  error  sources  considered  in  this  investigation. 
We  will  limit  the  discussion  here  to  multipath,  since  the  ionosphere  contribution 
will  be  examined  in  detail  in  Section  |222j  When  signals  from  a  GPS  satellite  reach  a 
receiver  via  multi (ple)-paths,  this  is  multipath.  Multipath  can  result  from  low  satellite 
elevation  and  poor  receiver  and  antenna  siting  with  respect  to  buildings,  terrain,  etc. 


As  noted  by  Hofmann-  Wellenhof  et  al.  2001  ,  in  severe  cases,  low-elevation  multipath 


can  result  in  loss  of  a  GPS  signal.  Less  severe  cases  can  result  in  measurement  errors. 
Mitigation  can  be  achieved  by  employing  choke-ring  type  antennas  with  better  gain, 
and  by  careful  siting  of  the  antenna.  Since  multipath  can  introduce  errors  and  outages, 
researchers  often  apply  a  mask  to  exclude  data  from  satellites  below  elevations  up  to 
40°  |  Thomas  et  al. ,  2004  ,  although  15°  seems  common  [see,  e.g.,  Datta-Barua  et  al. 


2003  .  For  our  research  we  applied  a  15°  mask  to  all  the  data. 

The  sum  of  the  measurement  errors  listed  above  is  called  the  User  Equivalent 
Range  Error  (UERE)  and  can  be  multiplied  by  a  factor  called  dilution  of  precision  to 
obtain  position  error. 


2.1.4  Dilution  of  Precision.  Dilution  of  Precision  (DOP)  directly  relates 
measurement  errors  to  position  errors.  To  arrive  at  DOP,  we  must  begin  with  the 


pseudo-range  equation  (2.2)  introduced  earlier.  Following  a  development  of  the  DOP 


equation  by  J.  Raquet  (EENG  533,  Navigation  Using  the  GPS,  course  handouts, 
2006),  and  also  Dempster  |2006],  recall  from  Subsection  2.1.1  that  four  pseudo¬ 
range  equations  are  necessary  to  solve  for  position.  These  equations  can  be  expanded 
using  a  Taylor  series,  linearized,  and  solved  through  iterative  techniques  involving  an 
approximate  solution  (i.e.  guess).  The  linearized  equations  can  be  expressed  as  a 
matrix: 


A p  =  HAx  or 

(2.3) 

Ax  =  H_1Ap 

(2.4) 

9 


where  Ax  is  the  vector  of  position  and  receiver  clock  corrections,  A p  is  a  vector 
of  range  error  (i.e.  difference  between  actual  range  and  guess),  and  H  is  the  4x4 
matrix  consisting  of  unit  vectors  (x,  y,  z)  between  the  satellite  and  the  guess  position, 
plus  a  clock-error  term.  When  more  than  four  satellites  are  used,  the  system  is 


overdetermined  and  can  be  solved  using  least-squares  techniques.  Equation  (2.4) 
becomes: 

Ax  =  (HtH)  _1  HrAp  (2.5) 


When  two  assumptions  are  made:  that  all  measurements  have  the  same  variance,  and 
the  measurement  errors  are  uncorrclated,  we  obtain  through  least-squares  theory: 


Cx  —  (HtH)  -1  (j2p 


(2.6) 


where  Cx  is  the  “covariance  matrix  of  calculated  position  and  clock  error”,  and  a2p 
is  the  variance  of  the  range  measurement.  If  the  H  matrix  is  transformed  into  local 


coordinates  (ENU),  the  first  term  on  the  r.li.s.  of  Equation  (2.6)  becomes: 


Di  i 

D 12 

D\3 

D 14 

D21 

D22 

D23 

D24 

D31 

D32 

D33 

D34 

-D41 

D42 

D43 

D44 

(2.7) 


where  Dij  are  in  the  local  (geodetic)  coordinate  frame,  indicated  by  the  G  (geode¬ 
tic)  superscript.  This  is  the  DOP  matrix,  and  it  directly  relates  position  errors  to 
measurement  errors  as  we  noted  in  the  beginning  of  this  section. 

DOP  can  be  thought  as  being  inversely  proportional  to  the  volume  of  a  body: 


.  .  .  This  body  is  formed  by  the  intersection  of  the  [four]  site-satellite  vec¬ 
tors  with  the  unit  sphere  centered  at  the  observing  site.  The  larger  the 
volume  of  this  body,  the  better  the  satellite  geometry.  Since  good  geome¬ 
try  should  mirror  a  low  DOP  value,  the  reciprocal  value  of  the  volume  of 
the  geometric  body  is  proportional  to  DOP.  The  critical  configuration  is 
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Hofmann- Wellenhof  et  al. 


given  when  the  body  degenerates  to  a  plane. 
2001 


DOP  variants  can  then  be  extracted  from  the  matrix.  Geometric  DOP,  or 


GDOP  is  the  square  root  of  the  trace  of  matrix  (2.7).  As  noted  by  J.  Raquet  (EENG 
533,  Navigation  Using  the  GPS,  course  handouts,  2006),  the  key  relationship  is 


GDOP  x  trUEKB  =  /<xf +  +  + 


(2.8) 


which  “relates  UERE  with  the  RSS  (root  square  sum)  of  the  errors.”  We  mentioned 
in  Chapter  [T]  that  our  research  is  concerned  with  HDOP,  the  Horizontal  Dilution  of 
Precision,  which  is  defined  as: 


HDOP  =  \/D  ii  +  D22.  (2.9) 

In  practice,  a  DOP  value  of  one  is  perfect,  10  is  generally  considered  fair-poor,  al¬ 
though  the  relative  badness  depends  on  the  particular  circumstances.  The  ionosphere 
influences  DOP  by  degrading  or  obliterating  GPS  transmissions,  thereby  reducing 
the  number  of  satellites  which  are  available  for  a  navigation  solution.  In  order  to 
understand  how  and  why  the  ionosphere  affects  GPS  transmissions,  we  need  to  know 
something  about  those  transmissions. 

2.1.5  Characteristics  of  GPS  Transmissions.  GPS  transmits  right-hand 
circularly  polarized  waves  on  two  carrier  frequencies,  designated  LI  and  L2,  at  1574.42 
MHz  and  1227.6  MHz  respectively.  The  corresponding  wavelengths  are  19.03  cm  for 
LI  and  24.42  cm  for  L2.  The  carriers  are  modulated  by  Pseudo-Random  Noise  (PRN) 
codes,  of  which  there  are  two  types,  Coarse- Acquisition  (C/A  Code)  code  and  Precise 
(P-Code)  code.  The  C/A  code  is  modulated  on  LI  only  and  is  available  to  civilian 
users,  while  the  P  code  appears  on  both  LI  and  L2,  and  although  unclassified,  is  not 
typically  transmitted.  Instead,  the  P  code  is  encrypted  with  a  classified  encryption 
key  and  called  Y  code,  or  P(Y)  code.  Each  satellite  is  assigned  its  own  PRN  code, 
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although  a  particular  PRN  can  be  reassigned  to  another  satellite.  The  PRN  codes  are 
identified  by  an  integer  in  the  range  of  1  to  32,  and  the  satellites  are  often  referred 
to  by  their  PRN.  In  addition  to  the  PRN  codes,  the  carrier  is  also  modulated  with  a 
navigation  message.  This  research  was  confined  to  reception  of  LI,  C/A  code  only. 
Much  can  happen  to  the  signal  once  it  leaves  the  satellite,  particularly  as  it  travels 
through  the  ionosphere,  the  subject  of  the  next  section. 


2. 2  The  Ionosphere 

Becoming  acquainted  with  the  ionosphere  is  the  next  step  toward  understanding 
scintillation  effects  on  GPS.  In  the  preceding  section,  the  Earth’s  ionosphere  was 
mentioned  as  a  source  of  error  in  GPS  location  estimates.  We  begin  this  section 
by  introducing  the  structure  and  characteristics  of  the  ionosphere.  A  description  of 
the  effect  of  Earth’s  magnetic  field  on  the  ionosphere  follows. We  conclude  with  a 
discussion  of  the  interaction  of  EM  waves  with  ionospheric  plasma. 


2.2.1  Characteristics  of  the  Ionosphere.  Earth’s  ionosphere  forms  in  re¬ 
sponse  to  the  interaction  of  energy  from  the  sun  with  neutral  atmospheric  constituents. 
Chemical,  electrodynamic,  kinematic  and  mechanical  processes  are  all  at  work  at  var¬ 
ious  latitudes  and  altitudes  and  within  the  ionosphere.  The  ionosphere  can  be  strati¬ 
fied  according  to  the  processes  (e.g.  photochemical,  transport)  products  (e.g.  species) 
or  phenomena  (e.g.  reflection  of  EM  energy,  density  peaks).  Characteristics  of  the 


ionosphere  are  shown  given  Table  2.1 


Table  2.1  shows  the  primary  processes  in  the  D,  E,  and  F'i  regions  are  pho¬ 
tochemical  (including  photo-ionization,  recombination,  charge  exchange,  etc),  while 
a  transition  to  transport  dominance  occurs  within  the  F2  layer.  The  Topside  iono¬ 
sphere  is  dominated  by  transport  mechanisms,  including  ambipolar  diffusion  as  well 
as  gravitational,  E  x  B,  and  neutral  wind  drifts.  Figure  |2.1  is  a  modeled  electron 
density  profile  with  ionospheric  layers  from  course  notes,  PH2514,  Naval  Postgradu¬ 
ate  School,  R.C.  Olsen,  2003;  details  are  included  in  the  caption.  Profiles  are  given 
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Tabic  2.1:  A  summary  of  mid  and  low-latitude  ionospheric  layers  and  character¬ 

istics  summarized  from  \S chunk  and  Nagy,  2000 1.  Note  the  transition  to  transport 


that  occurs  in  the  upper  F2  layer.  This  will  be  important  later  when  discussing  the 
movement  of  plasma  in  the  equatorial  ionosphere. 

Ionospheric  Layers  and  Characteristics 


Layer 

Height  km 

Dominant  Process 

Energy  A 

Dominant 

Species 

D 

60-100 

Photochemical 

UV  &  X-ray 

±  ions  and 

“water  clus¬ 
ter”  ions 

E 

100-150 

Photochemical 

EUV 

NO+  Oj  N+ 

Fi 

150-250 

Photochemical 

EUV 

0+ 

f2 

250-30CQ 

Photochemical- 

Transport 

EUV 

0+ 

Topside 

Iono¬ 

sphere 

>300 

Transport  (Diffusion, 

Drifts) 

EUV 

0+ 

“Height  of  Electron  Density  Peak 


for  both  day  (solid)  and  night  (dashed)  as  well  as  for  solar  maximum  (right  profile) 
and  solar  minimum  (left  profile).  Notice  how  the  layers  decay  to  varying  degrees  at 
night  in  the  absence  of  photo-ionization.  The  Fi  layer  erodes  entirely,  while  the  E 
and  F2  regions  remain,  sustained  by  transport  and  slow  recombination,  respectively. 
Additionally,  note  the  increase  in  electron  density  during  solar  maximum. 


2.2.2  Earth’s  Magnetic  Field  Interacts  with  Ionospheric  Plasma.  The  seeds 
of  scintillation  are  planted  in  the  E-region  of  the  ionosphere.  Here,  tidal  winds 


drive  an  eastward  dynamo  current  in  the  sunlit  hemisphere  Kelley  1989,  chap.  3]. 
The  resulting  eastward  electric  held  is  perpendicular  to  the  geomagnetic  held  lines, 
which  in  turn  are  virtually  horizontal  at  the  geomagnetic  equator.  Referring  to  Fig- 
2.2  de  La  Beaujardiere  et  al.,  2004]  (looking  east),  the  E  x  B  drift  drives  the 


ure 


plasma  upward.  The  plasma  then  diffuses  along  B  held  lines  under  the  influence 
of  gravity,  gn,  and  pressure  gradients,  Vyp^  as  shown  on  the  left  and  right  sides  of 


Figure  |2.2[  It  accumulates  ~15-20°  north  and  south  of  the  dip-equator.  This  phe¬ 
nomena  is  called  the  Equatorial  Fountain.  The  accumulation  of  plasma  poleward  of 
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Figure  2.1:  Electron  Density  Profile  with  Ionospheric  Layers  from  course  notes, 

PH2514,  Naval  Postgraduate  School,  R.C.  Olsen,  2003.  “The  International  Reference 
Ionosphere  1995  (IRI-95)  model  was  run  for  the  location  of  Monterey,  CA  (geographic 
latitude  =  36.5,  geographic  longitude  =  238),  for  July  4,  1989  (solar  maximum)  and 
July  4,  1995  (solar  minimum).  Calculations  were  done  at  OLT  (local  midnight),  and 
12  LT  (local  noon)”  Notice  how  the  Fi  layer  vanishes  almost  entirely  during  the  night. 
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the  geomagnetic  equator  is  known  as  the  Appleton  Anomaly  and  the  locations  are 


labeled  Ionization  Crests  in  Figure  2.2  A  meridional  neutral  wind,  Vn||,  that  blows 
from  the  winter  to  summer  hemisphere  causes  plasma  in  the  summer  hemisphere  to 
be  forced  upward.  Conversely,  a  downward  drift  occurs  in  the  winter  hemisphere  thus 


the  Appleton  Anomaly  is  asymmetric  across  the  equator.  Figure  222  depicts  neutral 
wind  but  not  the  resulting  asymmetry. 


The  processes  described  above  and  illustrated  in  Figure  |2.2|  are  confined  to 
the  daylight  hemisphere.  As  Earth  rotates  into  darkness,  dramatic  changes  take 
place  which  can  lead,  ultimately,  to  scintillation.  At  the  boundary  between  day  and 
night  (the  terminator),  the  eastward  electric  held  increases  in  response  to  increasing 
neutral  winds  [Schunk  and  Nagy ,  2000  and  “polarization  charges  within  conductivity 
gradients  at  the  terminator”  |  de  La  Beaujardiere  et  al. ,  2004]  before  becoming  westerly. 
The  upward  ExB  drift  intensifies,  raising  the  F  layer.  This  is  called  the  Pre-Reversal 
Enhancement  ( PRE).  As  darkness  falls,  photo-ionization  ceases  and  the  F  region  is 
quickly  eroded  from  below  as  recombination  takes  place.  The  combination  of  lift 
and  bottom-side  erosion  results  in  steep  upward  density  gradient  forces  at  the  base 
of  the  F  layer,  directly  opposed  by  gravity.  This  configuration  is  characterized  by 
Rayleigh- Taylor  Instability,  in  which  a  heavier  fluid  is  supported  by  a  lighter  fluid. 


If  perturbed,  the  fluid  turns  over,  as  illustrated  in  Figure  |2.3[  creating  low-density 
structures  (often  referred  to  as  bubbles)  that  rise  through  the  denser  fluid.  In  the 
equatorial  ionosphere,  they  are  known  as  Equatorial  Plasma  Bubbles  (EPBs)  although 
their  shape  more  closely  resembles  elongated  tubes. 


Figure  2.4  illustrates  perturbation  growth  in  two  dimensions.  The  F-region 
density  gradient  arising  from  the  pre-reversal  enhancement  of  the  eastward  electric 
held  is  modeled  by  a  step-function.  The  magnetic  held  is  directed  into  the  page  (look¬ 
ing  north),  gravity  is  directed  downward  and  the  density  gradient  is  opposite  gravity. 
A  sinusoidal  perturbation  is  imposed  on  this  arrangement.  Assuming  the  plasma  is 
virtually  collisionless,  gravitational  drift  results  in  a  current  J  —  nM g  x  B/B2  which 


is  directed  eastward  as  shown  on  the  left  of  Figure  2.4  Charges  accumulate  at  the 
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2004 


Figure  2.2:  Appleton  Anomaly  formation  after  de  La  Beaujardiere  et  al. 

The  Equatorial  Fountain  is  indicated  by  the  upward  E  x  B  drift  arrow  at  the  top  of 
the  figure.  Note  the  E  held  has  been  mapped  from  the  E-region  to  the  F-region  along 
the  magnetic  held  lines  which  act  as  equipotentials.  Plasma  is  then  driven  down  the 


held  lines  by  gravity  gy  and  pressure  gradient  V||p1 
wind,  Vn|j  is  also  shown. 


forces.  The  meridional  neutral 
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Figure  2.3:  Rayleigh  Taylor  Instability  image  from  http://web.arizona.edu/  fluid- 
lab /ray leigh.html.  A  dense  fluid  is  held  aloft  by  a  less  dense  fluid  below.  A  pertur¬ 
bation  can  trigger  overturning,  leading  to  the  production  of  lower  density  structures 
(bubbles)  that  rise  through  the  denser  fluid.  In  this  figure,  the  dense  fluid  is  white. 


inflection  point  of  the  sinusoidal  discontinuity  since  the  current  is  larger  (oc  n)  in  the 
dense  region.  In  response  to  the  charge  accumulation,  polarization  E  fields  develop 
as  shown.  Subsequent  E  x  B  drift  causes  the  less  dense  plasma  to  rise  and  the  dense 
plasma  to  sink,  amplifying  the  perturbation  and  eventually  leading  to  overturning 
and  bubble  formation. 


Using  linear  perturbation  theory  and  assuming  a  sinusoidal  perturbation  as 
above,  a  simple  expression  for  the  growth  rate  7 nr  of  the  EPBs  is  found  to  be: 


f  g_\  1  dnp 
\vin)  n0  dz 


(2.10) 


A  derivation  can  be  found  in  Kelley  1989,  chap:  4],  Gravity,  g  and  density 
gradient,  ///,  were  introduced  in  our  discussion  of  Figure  2.4,  although  here  we  are  no 
longer  considering  a  simple  step-function.  For  fastest  growth,  the  plasma  gradient 
should  be  positive  (dense  fluid  aloft)  and  large.  The  ion-neutral  collision  frequency, 


uin,  enters  through  the  Pederson  conductivity  (  Kelley  1989,  see])  and  should  be 
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Figure  2.4:  Perturbation  Growth  Kelley,  1989 


small,  as  it  is  at  high  altitudes.  Recall  the  pre-reversal  enhancement  drives  the  F¬ 


layer  upwards  near  sunset,  further  decreasing  vin.  Equation  (2.10)  accounts  for  the 


role  played  by  the  gradient  and  ion-neutral  collision  frequency,  but  other  forcing 
mechanisms  can  drive  the  instability,  including  an  eastward  electric  held  E  and  the 
neutral  wind  u„. 


Beginning  with  the  gravitationally  driven  drift  current  shown  in  Figure  [274],  and 
adding  terms  for  the  eastward  electric  held  and  the  neutral  wind,  we  obtain: 


nrriig  _  .  „ , 

J  ss  6  +  o,,E  +  ar  (u„  x  B) 


(2.11) 


Recalling 


<jp 


nrrii  uin 


B2 


and  substituting 


nrrii  q  vm  E  un  B 

J  &  — +  n  mi  ——5-  +  n  mi  vin - 


B 


B2 


B2 


(2.12) 
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when  the  ratio  of  cyclotron  frequency  to  species-neutral  collision  frequency,  k  for  ions 
Ki  1  and  Ke  Ki.  Further  simplifying, 


J 

J 


9  + 

g 

L'in 


Vin  E  l^in  Un  B 

B  B 
E  unB 

H - h  — — 

B  B 


(2.13) 

(2.14) 


Thus  the  first  term  of  (2.14)  is  in  a  form  equivalent  to  the  gravitational  term  in 


Equation  (2.10).  We  can  now  append  the  corresponding  electric  held  and  neutral 


wind  terms  to  Equation  (2.10): 


7  RT 

'JRT 


g  1  dn o  El  dn o  1  driQ 

Vin  no  dz  B  n0  dz  n  n0  dz 
1  9np  (_9_  Eo 

n0  dz  V  i',n  B  U 


(2.15) 

(2.16) 


An  eastward  (positive)  electric  held  and  a  downward  (positive)  neutral  wind  con¬ 
tribute  to  growth,  however  there  is  typically  no  strong  downward  neutral  wind  com¬ 
ponent.  The  neutral  wind  still  plays  a  role,  however,  if  the  plasma  layer  is  tilted 
with  respect  to  the  vertical.  This  can  happen  locally  in  response  to  the  pre-reversal 
enhancement.  The  addition  of  a  horizontal  component  of  Vn  allows  the  eastward 
neutral  wind  to  influence  the  stability.  Letting  a  refer  to  the  angle  the  gradient  Vn 
makes  with  the  horizontal  and  accounting  for  the  eastward  and  upward  components 
of  the  electric  held,  we  obtain: 


7-RT  — 


1  dn0 
no  dz 


—  cos  a  + 

V in 


rr'east 

-^0 

B 


COS  <7  + 


(£7  +  uenastB) 


B 


sin  er 


(2.17) 


Earlier  in  this  section  we  said  the  bubble  growth  occurs  as  the  Earth  rotates  into 
darkness.  During  the  day,  the  E-region  dynamo  overwhelms  the  F-region  dynamo, 
inhibiting  bubble  growth.  During  the  night,  however,  the  E-region  decays,  the  F- 
region  dynamo  dominates,  and  bubbles  can  grow.  To  capture  this  diurnal  behavior, 
we  include  the  weighted  average  of  the  field-line  integrated  Pederson  conductivities 
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Tjp  in  the  E  and  F  region: 


7  RT  — 


Yf 
u p 


1  drin 


Ep  +  Yf,  n0  dz 


g  Eenast  (E%p  +  uiastB) 

—  cos  a  H - cos  a  +  AA - — - -  sm  a 


B 


B 


(2.18) 


Equation  2.18  now  captures  the  diurnal  behavior  of  plasma  bubble  formation  and 
accounts  for  the  influence  of  the  density  gradient  AAhi  which  should  be  strong  and 
directed  upward,  the  eastward  electric  held  and  eastward  neutral  wind. 

EPBs  have  been  the  subject  of  research  for  decades,  although  they  were  mysteri¬ 
ous  phenomena  on  remote  sensing  imagery  before  the  cause  was  determined.  At  first, 
the  term  Equatorial  Spread-F  (ESF)  was  used  to  describe  the  ionosonde  returns  ob¬ 
served  when  EPBs  were  encountered.  The  reflected  echo  depicted  a  spread-F  region, 
either  in  frequency  or  range.  In  the  1970s,  returns  from  a  50  MHz  radar  at  Jicamarca, 
Peru  displayed  structures  referred  to  as  plumes.  In  the  subsequent  decades,  the  cause 
of  the  unusual  returns  was  identified  as  EPBs.  Much  of  what  has  been  learned  about 


EPBs  is  summarized  by  Schunk  and  Nagy  2000  and  more  recently  by  Dandekar  and 


Groves  2004  and  de  La  Beaujardiere  et  al.  2004  .  Table  2.2.2  presents  a  synopsis. 


The  number  of  EPBs  parallels  the  solar  cycle  and  peaks  at  the  equinoxes,  but  shows 
little  response  to  day  to  day  variations  in  solar  flux  or  geomagnetic  activity.  Seasonal 
activity  is  linked  to  the  Equinoxes  for  GPS  frequencies,  but  at  VHF  wavelengths  peak 
activity  differs  by  longitudinal  sector.  Individual  bubbles  begin  to  form  about  an  hour 
after  sunset  and  each  lasts  for  about  80  min.  The  bubbles  tend  to  form  every  160 
min,  with  peak  activity  around  2200L.  Two  to  three  bubbles  form  each  day.  The 
longitudinal  extent  of  the  activity  spans  1000s  of  kilometers.  The  bubbles  drift  east¬ 


ward  anywhere  from  100-200  m/s  Immel  et  al,  2003  ,  tilting  westward  with  height. 


Their  vertical  velocity  typically  ranges  from  100-500  m/s  although  velocities  can  be 
ten  times  faster.  They  tend  to  form  fossil-bubbles  (no  longer  actively  growing)  as  the 
night  wears  on  and  upward  drift  subsides.  The  plasma  density  can  drop  as  much  as 
102  from  the  outside  to  the  inside  of  an  EPB,  and  this  density  variation  is  critical  to 
scintillation. 
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Table  2.2:  Temporal  and  Spatial  Characteristics  of  Equatorial  Plasma  Bubbles 

&  Scintillation  based  on  remote  sensing  and  in-situ  data,  compiled  from 
and  Nagy  1 2000 


Schunk 


Dandekar  and  Groves  2004  ,  Immel  et  al.  2003  ,  and  K.  Groves 


(personal  communication,  2006). 


Temporal  Characteristics 

solar /seasonal  sensitivity 

parallels  solar  cycle  but  no  short  term 
(F10.7)  solar  flux  dependence;  worst 
at  Equinox, 

geomagnetic  sensitivity 

worst  under  Appleton  Anomaly  crests, 
less  at  geomagnetic  equator.  No  cor¬ 
relation  with  Kp  or  Dst 

individual  bubble  duration 

l:20h 

interval  between  formation 

2:40h 

#  bubbles  d-1 

2-3 

daily  cycle 

begins  ~  lh  after  sunset,  peaks  2200L, 
dissipates  by  sunrise 

Spatial  Characteristics 

longitudinal  extent  of  activity 

1000s  km 

longitudinal  EPB  interval 

10-1000  km 

horizontal  velocity 

100-200m/s 

vertical  velocity 

100-500  m/s,  40%  reach  500-5000  m/s 

vertical  extent 

«  100-1500  km 

trajectory 

drift  eastward,  tilt  westward  with 
height 

decay 

”  fossil”  bubbles  form  as  Spread- F  & 
upward  drift  ends 

density  variation 

up  to  102  lower  plasma  density  inside 
EPB 

21 


2.3  Scintillation:  A  Consequence  of  Ionospheric  Irregularities 

The  final  step  toward  understanding  scintillation  effects  on  GPS  culminates  with 
a  discussion  of  scintillation  itself.  We  will  define  scintillation  and  examine  its  influence 
on  EM  waves.  The  section  will  conclude  with  a  discussion  of  research,  measurement 
and  modeling.  First,  however,  we  must  understand  how  EM  waves  travel  through  a 
plasma. 

2.3.1  Wave  Propagation  in  Plasma.  For  simplicity,  consider  the  dispersion 
relation  for  EM  waves  propagating  in  a  cold,  unmagnetized,  isotropic  electron  plasma 
given  by: 


U 


2 


where  up  = 


2  I  2;  2 

up  +  c  k 

47 mee2 
me 


(2.19) 

(2.20) 


The  terms  are  the  frequency  of  the  transmitted  wave,  u,  the  plasma  frequency,  up, 
the  wave  number,  k,  and  the  number  density  of  electrons,  ne.  The  phase  velocity  can 
be  extracted  by  rearranging  the  terms: 


Vl  =  *2  =  #  + 


U 


oc 


(2.2!) 


Similarly,  the  group  velocity  is  given  by: 


2  -  duJ  I  1 

v  g  =  w  =  c[ 


T<c 


(2.22) 


where  kr  is  introduced  representing  the  real  portion  of  the  wave  number.  As  an  EM 
wave  passes  through  a  plasma,  the  phase  velocity  is  advanced,  and  the  group  velocity 
delayed.  This  results  in  overestimated  GPS  code  pseudoranges  and  underestimated 
carrier  pseudoranges.  Further  manipulations  involving  the  refractive  index  produce 
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the  following  relationships: 


K?°=^r  '"-d 


e  Oso 


/  nRd 


e  UsO 


(2.23) 


where  A  represents  the  difference  between  measured  and  geometric  range,  the  constant 
40.3  arises  from  a  Taylor  series  expansion  of  the  phase  refractive  index,  and  integration 
occurs  along  the  geometric  path  dso  | Hofmann- Wellenhof  et  al..  2001  .  The  quantity 


f  ne  dso  is  the  Total  Electron  ContentTEC.  The  important  point  from  Equation  (2.23 ) 
is  that  using  the  relationship  between  TEC  and  frequency,  the  ionospheric  refraction 
can  be  exploited  to  obtain  TEC  measurements  from  dual  frequency  receivers  and 
conversely,  the  same  relationship  can  be  modeled  to  remove  the  ionospheric  refraction 


from  single  channel  receivers  as  described  by  Klobuchar  Parkinson  et  al. ,  1996  . 


2.3.2  Definition,  Characteristics  and  Models.  The  situation  for  scintillation, 


sadly,  is  not  so  simple.  Groves  et  al.  1996  provide  a  concise  description  of  scintillation: 


Scintillation  of  radio  waves  occurs  when  the  waves  encounter  small-scale 
spatial  variations  in  refractive  index. . .  Because  the  scattering  and  refrac¬ 
tive  properties  of  the  medium  vary  along  an  incident  wave  front,  the  trans¬ 
mitted  wave  front  is  composed  of  a  mixture  of  varying  amplitude  and  phase 
and  an  irregular  interference  pattern  results. .  .  The  scintillation  observed 
on  the  ground  is  an  indication  of  the  rate  at  which  the  interference  pat¬ 
tern  drifts  past  the  observer,  determined  by  the  physical  drift  speed  of 
the  [EPBs],  the  velocity  at  which  the  source  wave  is  cutting  across  the 
ionosphere  due  to  satellite  motion,  and  the  motion  of  the  observer,  if  any. 


To  understand  the  relationship  between  refractive  index  J\f  and  plasma  density 


n,  we  return  to  the  dispersion  relation,  Equation  2.19 


2  2  i  2/2 

OJ  =  LVp  +  c  k 


The  refractive  index  is  defined  as  : 


a r2  = 


{u/kf 


(2.24) 
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Substituting  2.24  into  2.19  we  obtain 


A/-2  =  l-^ 

a r 


(2.25) 


Recalling  Equation  |2.20  and  substituting  in  2.25 


mu- 


Differentiating  with  respect  to  x  and  rearranging  terms: 


dA F  f  2  d  ne 

“4TO 


mu 


1  - 


47mPe2\  1//2 


mu- 


(2.26) 


(2.27) 


This  shows  that  the  change  in  refractive  index  over  a  distance  is  directly  proportional 


to  the  change  in  electron  density  over  that  same  distance.  Recall  from  Table  2.2.2 


that  the  plasma  density  can  change  by  two  orders  of  magnitude  across  an  EPB. 
Figure  [275]  is  an  artistic  interpretation  of  the  interference  pattern  caused  by  the  EPBs. 
As  EM  waves  pass  through  ionospheric  density  depletions,  the  waves  are  refracted  and 
scattered,  leading  to  “patches”  of  focused  and  defocused  energy.  The  arrow  shows 
these  “patches”  drifting  past  the  receiver  (person)  located  on  the  ground.  The  final 
result  is  fluctuations  in  the  intensity  of  the  received  signal 

Researchers  have  historically  modeled  scintillation  using  a  one-dimensional  phase 
screen.  A  1-D  phase  screen  can  be  thought  of  as  an  arrangement  of  lenses  alternately 
focusing  or  defocusing  the  incident  EM  energy  as  it  passes  through  the  screen.  The 
screen  moves  relative  to  the  transmitter  and  receiver  as  described  above.  A  simple 
model  of  of  a  discrete  phase  screen  is  described  by  Knight  |2000|  depicted  in  Figure [2Tj 


In  the  figure,  the  shaded  circles  labeled  /  represent  “rod  like  lenses  aligned  to 
the  Earth’s  magnetic  field”,  dx  describes  a  portion  of  the  phase  screen  and  r  is  the 


distance  to  the  elements  Knight  2000,  Appendix  A],  He  gives  the  complex  amplitude 
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FROM  SATELLITE 


Figure  2.5:  Illustration  of  scintillation.  As  EM  waves  pass  through  ionospheric 

density  depletions,  the  waves  are  refracted  and  scattered,  leading  to  “patches”  of 
focused  and  defocused  energy.  The  arrow  shows  these  “patches”  drifting  past  the 
receiver  (person)  located  on  the  ground.  The  final  result  is  fluctuations  in  the  intensity 
of  the  received  signal. 


Incident  wave 


Figure  2.6: 

the  shaded  circles  labeled 


2000 


Illustration  of  a  thin  Phase-Screen  model  of  scintillation  by  Knight 

rod  like  lenses  aligned  to  the  Earth’s 


magnetic  field” 
the  elements. 


I  represent 

da;  describes  a  portion  of  the  phase  screen  and  r  is  the  distance  to 
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of  the  signal  as: 


I  = 


—j=  ( 1  —  2  /  sin(P0  +  ^>/2)  sin$/2  dx  ) 
Vr A  V  Jir  ) 


g  =  Ax(2. 


cos 


(P0  +  $/2)  sin <f>/2  da: 


(2.28) 

(2.29) 


where  /  and  Q  represent  the  in-phase  and  quadrature  components  of  the  signal,  $  is 
the  phase  perturbation  produced  by  the  lenses,  and  P0  =  — 7t/4  —  27r(r  —  h)/  A,  and 
the  quantity  \fr\  is  the  Fresnel  radius,  zj.  The  Fresnel  radius  indicates  whether  an 
ionospheric  irregularity  will  have  an  impact  on  amplitude.  If  the  scale  of  the  irreg¬ 
ularity  is  much  larger,  amplitude  variations  are  minimal.  At  or  below  Zf,  amplitude 
variations  are  significant.  Variations  are  also  significant  if  strong  density  gradients 


are  present  Knight,  2000,  Appendix  A.].  Beach  et  al.  2002  studied  the  phase  screen 
model  and  found  it  generally  able  to  produce  an  estimate  of  the  magnitude  of  ampli¬ 
tude  scintillation  within  90%  of  its  true  value  when  “only  a  few  (~5)  evenly  spaced 
samples  per  Fresnel  radius”  are  obtained.  Zf  for  GPS  frequencies  is  «  300-400  m. 
The  severity  of  amplitude  scintillation  is  expressed  in  terms  of  the  S4  index,  where 


S4  = 


V(i2)  -  (I)2 
(I) 


(2.30) 


In  Equation  (2.30)  /  is  the  signal  intensity  measured  at  the  receiver  \Fremouw  et  al. 


1980  ,  or  the  square  of  the  amplitude.  Equation  2.30  is  simply  the  signal  intensity 


standard  deviation  divided  by  its  mean,  with  sampling  rates  and  averaging  periods 
appropriate  for  the  signal  being  characterized,  and  is  valid  for  weak  to  moderate  levels 


of  scintillation.  Davies  1990  presents  a  model  of  the  S4  index  based  on  physical 
parameters: 


S4  =  B  Xs  Z  L  K  sec2  y  (An2)  (2.31) 

where  B  is  a  parameter  that  depends  on  “geometrical  factors,  fundamental  constants, 
etc”,  A  is  the  wavelength  of  the  carrier,  Z  and  L  are  the  layer  thickness  and  height, 
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respectively,  x  is  the  angle  between  the  vertical  and  the  ray,  and  (Arig)  is  the  mean 
square  deviation  of  the  electron  number  density.  As  well  as  directly  linking  S4  to 
electron  number  density,  this  equation  shows  S4  increasing  as  the  path  nears  the 
horizon. 

Thus  far,  we  have  linked  EPB-induced  plasma  density  variations  to  amplitude 
fluctuations  that  can  be  characterized  by  the  £4  index.  This  development  can  be 
taken  a  step  further  by  relating  the  £4  index  to  the  amount  of  signal  fade  at  the 
receiver.  This  is  accomplished  through  the  Nakagami-m  probability  density  function: 

P(A)  d  A  =  d^r(ro)(^2)m  exp  (-m^  )  (2.32) 

where  T  is  the  gamma-function,  A  is  the  amplitude,  and 


m—  1  /  S4 


(2.33) 


Wheelon  2003,  chap.  10]  points  out  that  “m  =  1  corresponds  to  a  Rayleigh  distribu¬ 
tion  and  rn  —  1/2  corresponds  to  a  one-sided  Gaussian  distribution.”  He  notes  the 
fit  of  S4  and  signal  fade  to  the  Nakagami-m  distribution  is  empirical  and  not  derived 


from  Erst  principles.  Figure  2.7  depicts  the  Nakagami-m  PDF  where  “The  curves 
are  each  normalized  so  that  the  maximum  likelihood  corresponds  to  unit  probabil¬ 


ity”  j  Wheelon ,  2003  .  The  take-away  is  as  m  becomes  smaller  (and  hence  S4  becomes 


larger)  the  range  of  likely  amplitude  values  increases.  Finally,  Knight  2000  notes  the 


£4  index  scales  as  £4  oc  f(~p+3)/4  where  /  is  the  carrier  frequency,  p  is  the  spectral 
index  (typically  2.5  at  equatorial  latitudes),  when  S4  <0.5  (weak  to  moderate  scin¬ 
tillation).  By  the  scaling  equation,  the  GPS  L2  frequency  (which  is  lower)  will  be 
more  susceptible  to  scintillation  than  the  LI  frequency,  by  a  factor  of  1.4  according 
to  Knight  2000  .  This  has  implications  for  dual  frequency  GPS  users:  if  scintillation 


is  strong  enough  to  cause  outages  on  LI,  L2  will  probably  be  unavailable  as  well. 
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< 

CL 


8.686  log,,  (A/A^  )  in  dB 


Wheelon 


2003 


Figure  2.7:  Nakagaim-m  Probability  Density  Function  by 

amplitude  and  ARM$  is  the  root- mean-square  amplitude.  “The  curves  are  each  nor 


A  is 


malized  so  that  the  maximum  likelihood  corresponds  to  unit  probability ”  Wheelon 


2003  .  Notice  as  m  becomes  smaller  (and  hence  S4  becomes  larger)  the  range  of  likely 


amplitude  values  increases. 
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Thus  far,  we  have  reviewed  the  characteristics  of  GPS  and  the  ionosphere. 
Within  the  ionosphere,  we  have  examined  the  development  of  EPBs  and  their  role  in 
scintillation.  We  have  seen  how  the  changes  in  refractive  index  arising  from  plasma 
density  gradients  within  the  EPBs  alter  EM  waves  passing  through  them,  and  thus 
alter  the  GPS  signals.  We  saw  that  the  amplitude  fluctuations  caused  by  scintillation 
can  be  characterized  by  the  S4  index.  Finally,  we  were  able  to  relate  the  S4  values 
to  levels  of  signal  fade  at  a  receiver  through  the  Nakagami-m  probability  density 
function.  The  last  portion  of  Chapter  [XT]  will  be  devoted  to  a  review  of  some  of  the 
literature  published  about  scintillation  and  GPS. 


2.3.3  Research  on  Scintillation  and  its  Impacts  on  GPS.  Literature  on 
scintillation  and  GPS  can  be  divided  into  observational  and  modeling  studies  and 
comparisons  of  the  two.  In  the  past  decade,  a  network  of  GPS  receivers  to  monitor 
scintillation,  as  an  enhancement  to  the  Scintillation  Decision  Aid  (SCINDA)  |  Groves 


et  al. ,  1997  ,  was  deployed  throughout  the  equatorial  ionosphere.  Consequently,  large 


data  sets  have  become  available  for  analysis.  Additionally,  in-situ  measurements  of 
plasma  irregularities  have  continued  and  will  expand  with  the  anticipated  launch  of 
the  Communication/Navigation  Outage  Forecasting  System  (C/NOFS)  [ de  La  Beau 
jardiere  et  al,  2004  . 


Groves  et  al.  1996  presented  some  hirelings  regarding  scintillation  and  L-band 


satellites  through  1996.  Early  data  from  GPS  receivers  indicated  signal  fades  exceed¬ 
ing  15  dB  were  possible  and  more  than  one  GPS  satellite  could  potentially  be  dropped 
during  a  scintillation  event.  His  follow-up  1997  paper  |  Groves  et  al .,  1997  described 
the  SCINDA  network  and  the  temporal  and  spatial  variations  of  L-band  scintillation 
in  the  Atlantic  and  Americas  sector,  along  the  magnetic  dip  equator  as  well  as  under 


the  anomaly  crest.  This  information  has  already  been  presented  in  Table  |2.2.2|  and 
will  not  be  repeated  here. 

By  2001,  the  SCINDA  network  had  grown  and  with  it,  the  knowledge  of  scintil¬ 
lation  impacts  on  GPS.  Additionally,  the  solar  cycle  had  reached  its  peak,  providing 


29 


additional  fodder  for  comparisons.  Thomas  et  al.  2001  presented  findings  from  the 


Pacific  sector.  In  addition  to  annual  summary  plots  of  scintillation  levels  at  Marak 
Parak  and  Parepare  showing  the  seasonal  and  diurnal  peaks,  he  presented  regional 
correlation  data  that  showed  if  any  satellite  in  the  region  experienced  £4  >0.6  on  at 
least  one  link  during  the  night,  “there  was  a  mean  probability  of  75%  that  the  other 
stations  in  the  network  will  have  at  least  one  link  which  also  records  scintillation 
above  £4  =0.6  on  that  night”.  He  also  presented  the  number  of  satellite  links  sur¬ 
viving  during  scintillation  episodes,  finding  that  up  to  40%  of  the  links  were  affected 
under  moderate  scintillation  (£4  >0.6),  and  over  90%  for  weak  scintillation  (£4  >0.3), 
markedly  higher  than  earlier  estimates.  In  his  estimation,  this  not  severe  enough  to 
significantly  hamper  navigation  because  the  episodes  were  short  lived. 


Datta-Barua  et  al.  2003  compared  the  scintillation  impacts  on  single  and 


dual-frequency  GPS  users.  They  found  both  susceptible  to  impacts  from  equatorial 
scintillation,  with  the  L2  frequency  more  prone  to  interruptions.  They  also  note  if  L2 
is  lost,  dual  frequency  users  can  substitute  the  ionospheric  model  available  to  single 
frequency  users  and  essentially  obtain  the  same  precision  and  availability  as  LI  users. 


Thomas  et  al.  2004  studied  GPS  precision  degradation  in  the  Pacific  sector  caused 


by  scintillation  during  the  1998-2002  portion  of  the  solar  cycle.  Errors  up  to  30  m 
were  discovered,  with  an  “underlying  seasonal  modulation  from  about  5  m  to  15  m” . 

At  the  same  time  SCINDA  was  gaining  momentum,  serious  efforts  at  mod¬ 
elling  ionospheric  scintillation  produced  the  WideBand  MODel  (WBMOD)  in  1995. 


WBMOD  is  based  on  climatology  and  the  one-degree  phase  screen  model  Rino 


1979b  Naturally,  comparisons  between  predicted  and  observed  conditions  began. 


In  1999,  Knight  et  al.  1999  compared  predictions  with  observations  in  the  Pacific 


sector.  They  found  good  agreement  between  modeled  and  predicted  observations. 


A  similar  study  by  Cervera  et  al.  2001  supported  the  earlier  work  for  periods  of 


low  sunspot  activity,  except  WBMOD  ended  scintillation  prematurely  at  night.  They 


XA  description  of  this  current  model  can  be  found  at  http://www.nwra- 
az.com/ionoscint/wbmod.html. 
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found  WBMOD  performed  poorly  during  solar  maximum  at  the  equatorial  stations, 
underestimating  scintillation  in  some  cases  by  an  order  of  magnitude.  WBMOD  also 
predicted  the  end  of  nightly  scintillation  two  hours  earlier  than  was  observed.  Since 
many  weapons  systems  rely  on  GPS,  a  logical  step  was  to  investigate  the  impact  of 
scintillation  on  various  weapons  systems.  In  a  modeling  study  of  B-1B  and  F-15E 
weapons  platforms,  \Evans  2005  demonstrated  that  scintillation  degrades  accuracy 
and  can,  if  occurring  in  conjunction  with  jamming,  render  GPS  unavailable. 


In  the  intervening  years  since  Thomas  et  al.  2001  ,  more  has  been  learned  about 
scintillation  and  GPS.  Thomas’  original  exploration  and  particular  presentation,  how¬ 


ever,  was  never  expanded  to  include  other  sites  or  years.  Thomas  et  al.  2001  wrote 
“It  is  desirable  to  have  a  more  rigorous  statistical  analysis  of  the  S4  database  for 
specific  S4  thresholds,  for  example,  as  a  function  of  local  solar  time  and  day  or  month 
of  the  year.”  Our  work  takes  up  this  challenge  for  additional  years  and  longitudi¬ 
nal  sectors.  Additionally,  we’ll  examine  the  influence  of  imputation  on  the  resulting 
climatology,  and  take  the  analysis  a  step  further  by  calculating  HDOP  based  on  the  se¬ 
quential  removal  of  satellites  experiencing  given  levels  of  scintillation.  The  remaining 
chapters  are  devoted  to  this  analysis. 
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III.  Method 


Creating  and  analyzing  statistics  first  requires  data.  We  begin  this  chapter  by 
describing  the  source  and  characteristics  of  the  data  we  used.  Data  must  also 
be  scrutinized  to  remove  undesirable  artifacts.  We  outline  the  preparation  of  the  data 
and  discuss  the  events  surrounding  missing  data.  Finally,  we  detail  the  processing  of 
the  data. 


3.1  Data  Description 

We  examined  SCINDA  data  provided  by  the  Air  Force  Research  Laboratory 
(AFRL).  Keith  Groves  (2006,  personal  communication)  suggested  beginning  our  anal¬ 
ysis  with  data  from  Ascension  during  the  last  solar  maximum  (2001).  Its  position  un¬ 
der  the  Appleton  Anomaly  crest  should  provide  a  rich  scintillation  environment.  The 
Pacific  sector  station  used  for  comparison  was  Parepare  during  2000,  also  under  the 
anomaly  crest.  We  also  analyzed  data  from  Ancon,  which  lies  close  to  the  magnetic 
dip  equator  and  a  corresponding  Pacific  sector  station,  Marak  Parak,  both  from  2001. 
Data  from  the  off  solar  peak  years  1999  and  2002  at  Ancon,  Ascension,  Marak  Parak, 
and  Parepare,  as  well  as  1998  at  Marak  Parak  and  Parepare,  and  2003  at  Ancon  and 
Ascension  were  also  examined.  Characteristics  are  shown  in  Table  13.11  Site  locations 


are  shown  in  Figure  3.1  (a)  and  (b). 


Table  3.1:  SCINDA  Site  Characteristics.  Negative  signs  indicate  South  Latitude 
or  West  Longitude.  Magnetic  data  from  International  Geomagnetic  Reference  Field 
(IGRF10)  and  are  approximately  centered  on  the  data  period. 


Location 

ID 

Lat  ° 

Lon  ° 

Years 

Elev  m 

Dip  Z  ° 

Ancon 

GB 

-11.78 

-77.15 

98-04 

53 

1.57 

Ascension 

GE 

-7.94 

-14.39 

99-04 

160 

-39.87 

Marak  Parak 

GV 

6.31 

117.40 

97-02 

34 

-3.13 

Parepare 

GW 

-3.98 

119.65 

97-03 

65 

-25.63 
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Figure  3.1:  Scintillation  Network  Decision  Aid  (SCINDA)  sites  in  the  (a)  Atlantic 
&  Americas  and  (b)  Pacific  sectors  used  in  this  research.  The  operational  years  are 
included  in  the  lower  left  quadrant  of  each  location 
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Two  types  of  GPS  receivers  were  used  to  collect  data  at  these  locations,  a 


NovAtel  11-channel  narrow  correlator  receiver  described  by  Dierendonck  et  al.  1996 


or  an  Ashtech  ZY-12  described  by  Milner  2002  .  Signals  were  sampled  at  50  Hz  and 
20  Hz  respectively  and  recorded  each  minute.  An  example  of  the  raw  data  record  is 
contained  in  an  appendix.  The  following  data  were  recorded: 

•  Receiver  Latitude  and  Longitude  (Decimal  °) 

•  Year 


•  Day  of  Year 

•  Time  (s) 

•  PRN 

•  S4 

•  SV  elevation  &  Azimuth  (Decimal  °) 

•  SV  Latitude  &  Longitude  (Decimal  °) 

•  Assumed  300  km  Ionospheric  Pierce  Point  Latitude  &  Lon¬ 
gitude 


3.2  Data  Preparation  and  Irregularities 

FORTRAN  code  written  by  C.  Leakeas  (personal  communication,  2006)  was 
used  to  reformat  the  data  from  its  original  tabular  form  to  a  columnar  form  more 
amenable  to  Matlab®  ingest.  During  reformatting,  all  records  from  satellites  below 
15°  were  excluded  as  a  first  attempt  at  mitigating  the  effects  of  multipatbQ  The 
remaining  observations  presented  two  irregularities.  The  first,  and  easiest  to  reconcile, 
was  an  S4  value  of  9.00  that  appeared  periodically.  This  was  identified  as  “bad  data” 
in  the  original  record.  According  to  Ron  Caton  (personal  communication,  2007)  9.00 
was  ascribed  during  extreme  levels  of  scintillation.  Under  extreme  scintillation,  early 
receivers  would  produce  nonphysical  S4  values  and  these  were  replaced  with  9.00. 


1see  Subsection 


2.1.3 
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When  the  elevation  of  the  particular  satellite  was  less  than  zero,  £4  of  9.00  was 
ascribed.  This  happened  when  incorrect  ephemeris  data  was  linked  to  the  £4  data 
months  after  the  £4  data  had  been  collected.  Since  we  applied  a  15°  mask  in  an 
effort  to  remove  multipath,  the  remaining  high  values  could  be  attributed  to  high 
levels  of  scintillation.  Typically,  high  £4  values  surrounded  the  occurrence  of  9.00, 
again  supporting  this  conclusion.  The  receivers  also  recorded  1  <  £4  <  9.00  during 
extreme  levels  of  scintillation.  During  particularly  bad  scintillation  conditions,  drop¬ 
outs  would  occur  and  a  particular  satellite’s  £4  value  would  not  be  recorded.  This 
missing  data  presented  another  challenge  that  until  now  had  not  been  addressed. 


Data  gaps  occurred  on  time  scales  from  minutes  to  months.  It  could  affect  as 
few  as  one  to  as  many  as  all  available  satellites  reporting  in  a  given  minute.  Be¬ 
fore  attempting  to  understand  the  data  that  were  present,  it  was  necessary  to  study 
the  circumstances  surrounding  the  data  that  were  absent.  We  originally  limited  our 
analysis  to  a  time  scale  of  minutes  and  to  those  missing  records  bounded  by  complete 
records.  We  began  by  determining  how  much  data  were  missing. 


In  the  introduction  to  his  text  on  analyzing  missing  data,  Schafer  [l997|  says 
“When  the  incomplete  cases  comprise  only  a  small  fraction  of  all  cases  (say,  five 
percent  or  less)  then  case  deletion  may  be  a  perfectly  reasonable  solution  to  the 
missing-data  problem.”  This  assumes  the  data  are  missing  at  random.  Examining 
2001  data  for  Ancon,  Ascension,  and  Marak  Parak,  the  fraction  of  bounded  data 
missing  was  three  percent  or  less  for  all  but  Ascension,  which  was  approximately 
eight  percent.  Ignoring  the  data  seemed  plausible  at  three  of  four  sites.  However,  a 
significant  fraction  of  the  missing  data  might  be  associated  with  strong  scintillation, 
in  other  words,  Missing  Not  At  Random  (MNAR)[^]  Ignoring  the  those  data  could 
underestimate  scintillation  impacts. 


2For  definitions  of  missing  data  mechanisms,  and  a  good  introduction  to  deal¬ 
ing  with  missing  data,  see  James  Carpenter  and  Mike  Kenward’s  website  at 
http://www.lshtm.ac.uk/msu/missingdata/index.html 
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We  studied  characteristics  of  the  outages  in  2000  and  2001  in  an  attempt  to  sep¬ 
arate  those  related  to  scintillation  from  those  related  to  other  causes  (e.g.  multipath, 


equipment  failure,  etc.).  Figure  3.2  shows  the  cumulative  distribution  of  bounded 
data  outages  in  2001  for  Ascension  and  Ancon.  In  all  cases,  85%  or  more  of  the 
outages  lasted  10  min  or  less.  This  suggests  most  outages  are  associated  with  brief 
interruptions  associated  with  scintillation. 

We  also  examined  the  relationship  of  outage  duration  to  the  day  of  year.  Fig¬ 
ure  3.3  shows  this  relationship  for  Ascension  (a),  Ancon  (b),  Parepare  (c)  and  Marak 


Parak  (d).  The  number  of  long-duration  outages  peaks  at  Ascension  near  the  equinoxes, 


indicating  a  link  to  scintillation  (recall  Table  2.2.2).  This  behavior  was  less  evident  at 


the  other  stations.  The  longer  duration  outages  exhibited  some  structure  across  days. 
They  were  particularly  prominent  at  Parepare  (c),  but  can  also  be  seen  at  Marak 
Parak  (d)  and  Ascension  (a).  We  suspected  satellites  were  skirting  the  15°  masl^on 
these  days,  leading  to  long-term  outages. 

To  confirm  our  suspicions,  we  examined  a  plot  of  elevation  angle  just  prior 
to  an  outage  versus  outage  duration  for  these  stations.  The  results  are  shown  in 


Figure  3.4  For  the  Pacific  longitudinal  stations  at  Parepare  (c)and  Marak  Parak  (d), 
the  evidence  is  clear  that  the  longer  duration  outages  were  associated  with  elevation 
angles  below  20°.  Additionally,  the  gentle  rise  in  outage  duration  from  high  elevation 
angles  to  elevation  angles  just  above  20°  hints  that  scintillation  susceptibility  increases 
as  the  satellite  nears  the  horizon.  This  is  consistent  with  Equation |2.31| which  showed 
S*4  increasing  as  the  x  angle  approached  90°.  This  is  also  intuitively  satisfying  since 
the  possibility  of  multipath  increases  as  the  path  nears  the  horizon.  If  multipath  is 
excluded,  the  path-length  through  the  ionosphere  is  greater  at  low  elevations  (high  x) 
and  increases  the  probability  of  transiting  a  depleted  region.  Although  the  correlation 
for  elevations  less  than  20°  at  Ascension  (a)  is  less  striking,  the  general  trend  of  longer 
outages  associated  with  lower  elevation  angles  appears  again. 


3See  Subsection  2.1.3 
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Figure  3.2:  Distribution  of  Outage  Durations  at  Parepare  and  Marak  Parak  in 

2000,  and  Ascension  and  Ancon  in  2001.  Most  outages  ten  minutes  or  less  suggesting 
scintillation  as  the  culprit. 


Finally,  Figure  |T5| is  a  plot  of  satellite  elevation  (black  lines)  and  S4  (blue  lines) 
for  a  particular  GPS  satellite  (PRN  28)  for  day  126  at  Ascension  in  2001  (Ron  Caton, 
personal  communication,  2007).  The  abscissa  in  this  plot  represents  time  centered  on 
approximately  2200  UT,  the  ordinate  on  the  left  is  S4  level,  on  the  right  is  elevation  in 
degrees.  Note  the  dotted  line  at  the  bottom  of  the  figure  represents  the  10°  elevation 
mask.  I11  the  center  of  the  figure,  in  conjunction  with  the  spike  in  S4  values,  note  the 
elevation  has  slowly  decreased  to  near  the  10°  mask.  Since  we  were  applying  a  15° 
mask  to  our  data,  this  would  appear  as  an  extended  data  gap. 


Inspection  of  the  distribution  of  outages  for  a  given  elevation  angle  before  and 
after  an  outage  and  whether  the  satellite  was  ascending  or  descending  produced  the 


results  shown  in  Figure  3.6.  In  the  Atlantic  &  Americas  sector  at  Ascension  and  An¬ 
con,  the  link  between  low  elevation  angles  and  outages  is  readily  apparent.  The  link 
is  less  obvious  in  the  Pacific  sector  at  Parepare  and  Marak  Parak,  but  still  present. 
Another  interesting  feature  present  at  all  sites  but  Ancon  is  a  cluster  of  outages  that 
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Figure  3.4:  Elevation  angle  versus  outage  duration  for  Ascension  (a)  and  Ancon 

(b)  in  2001,  and  Parepare  (c),  and  Marak  Parak  (d)  during  2000  with  15°  elevation 
mask  applied. 


Figure  3.5:  Satellite  Elevation  Near  Mask.  Plot  of  satellite  elevation  (black  lines) 
and  S4  (blue  lines)  for  a  particular  GPS  satellite  (PRN  28)  for  day  126  at  Ascension 
in  2001.  The  abscissa  in  this  plot  represents  time  in  half-hour  increments  centered  on 
approximately  2200  UT,  the  ordinate  on  the  left  is  £4  level,  on  the  right  is  elevation  in 
degrees.  Note  the  dotted  line  at  the  bottom  of  the  figure  represents  the  10°  elevation 
mask.  (Ron  Caton,  personal  communication,  2007) 
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began  between  40°  and  60°  as  the  satellite  is  ascending.  Comparison  with  outage 


duration  in  Figure  |3.4|  suggests  these  high  elevation  angle  outages  are  short  lived, 
and  thus  more  likely  to  be  associated  with  scintillation.  The  long-duration  low  ele¬ 
vation  outages  are  more  likely  a  consequence  of  satellite  geometry  (e.g.  multipath, 
obstructions)  and  mask  angle. 

At  this  point  we  decided  to  examine  the  effect  of  excluding  the  outages  below  20° 
elevation.  By  doing  so,  we  were  able  to  eliminate  the  structures  associated  with  long 


duration  outages,  as  seen  in  Figure  |3.7[  Although  the  structures  had  vanished,  some 
long  duration  outages  remained.  We  looked  at  outages  lasting  over  10  min  at  Marak 
Parak  (Figure  [3~T  (d))  in  2000  that  persisted  even  after  filtering  the  satellites  below 


20°.  The  characteristics  of  those  outages  are  shown  in  Table  3.2  The  first  column  of 
the  table  indicates  the  UTC  day  of  year  in  which  the  outage  occurred.  The  duration 
of  the  outage  in  seconds  is  shown  in  the  second  column.  Columns  3-5  indicate  the 
maximum  S 4  before,  during  and  after  the  outage  as  reported  by  the  available  GPS 
satellites.  The  final  two  columns  indicate  the  elevation  angle  of  the  affected  satellite 
immediately  before  and  after  the  outage.  The  first  five  outages  all  occurred  during 
benign  scintillation  conditions.  Additionally,  the  outage  durations  are  identical,  the 
elevation  angles  are  nearly  so,  and  the  outages  occurred  on  five  consecutive  days. 
These  facts  point  towards  an  outage  mechanism  other  than  scintillation.  Of  the 


remaining  long  duration  outages  in  the  Table  T2  the  one  on  day  50  also  began  during 
quiet  scintillation  conditions,  with  a  high  S4  in  the  midst  of  the  outage,  followed  by 
low  S4  at  the  end.  The  low  £4  values  bordering  the  outage  cast  doubt  on  scintillation 
as  the  culprit,  despite  the  anomalously  high  £4  value  during  the  outage.  Outages  on 
days  89,  110,  123  and  292  all  began  near  the  20°  threshold  or  ended  near  or  below  it, 
raising  the  specter  of  multipath  or  mask  problems.  Thus  we  are  left  with  days  88  and 
285  as  potential  scintillation  candidates.  These  two  days  are  near  the  equinox  and  so 
scintillation  seems  more  likely,  although  the  day  88  outage  begins  with  a  maximum 
reported  £4  of  0.2,  again  indicating  quiescent  conditions  at  onset.  So  of  these  12  long 
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Frequency  Frequency 


Elevation  Angle  (-  indicates  descent)  Elevation  Angle  (-  indicates  descent) 


(a) 


(b) 


Elevation  Angle  (-  indicates  descent)  Elevation  Angle  (-  indicates  descent) 


(c) 


(d) 


Figure  3.6:  Elevation  Angle  Distribution  Before  and  After  Outage  for  Ascension 

(a)  and  Ancon  (b)  in  2001,  and  Parepare  (c),  and  Marak  Parak  (d)  during  2000  with 
15°  elevation  mask  applied.  Negative  values  indicate  descent. 
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duration  outages,  2  can  be  reasonably  attributed  to  scintillation  to  the  exclusion  of 
other  causes. 


Finally,  we  looked  at  the  correlation  between  cr|4  and  S±max.  Figures [3~8|  (a,  c,  e) 
show  the  correlation  one  minute  before,  during,  and  one  minute  after  outages  during 


contains  the  corresponding  correlation  coefficients.  The  high  degree  of  correlation 
between  the  maximum  S4  and  S4  variance  before  and  after  the  outages  suggests  that 
GPS  constellation  geometry  and  the  scale  of  EPBs  result  in  only  a  few  satellites  being 
scintillated  at  once.  In  other  words,  scintillation  events  typically  affect  one  or  two 
satellites  during  a  given  minute,  leading  to  a  large  S4  variance.  The  relationship 
between  the  maximum  S4  and  S4  variance  described  above  may  be  useful  to  impute 
missing  S4  values  within  a  given  minute,  i.e.  choosing  a  replacement  S4  value  such 
that  the  correlation  between  the  maximum  S4  and  cr|4  is  maintained.  Similar  statistics 
have  been  used  to  distinguish  between  multipath  and  scintillation  \Datta-Barua  et  al. , 
2003  . 


2001  at  Ancon,  while  Figures  3.8  (b,  d,  f)  show  the  same  for  Ascension.  Table  3.3 


In  addition  to  the  single-satellite  outages,  outages  involving  all  satellites  occa¬ 
sionally  occurred.  Simultaneous  outages  accounted  for  up  to  approximately  34%  of 


all  outages  as  shown  in  Table  3.4  (a)  at  Parepare.  However,  S4  was  less  than  0.3  in 
85-90%  of  the  simultaneous  outages,  as  shown  in  Table  3H]  (b).  Thus  the  majority 
of  the  simultaneous  outages  can  be  attributed  to  phenomena  other  than  scintillation. 
Interestingly,  however,  the  duration  of  simultaneous  events  was  almost  always  two 
minutes  or  less. 


Returning  to  the  present  research,  although  only  four  sites  were  sampled  we  are 
confident  the  following  characteristics  are  applicable  in  general: 


•  Most  outages  are  short,  lasting  less  than  ten  minutes  (Figure  [372]). 

•  Outages  durations  reflect  to  some  degree  the  seasonal  peaks  in  scintillation 
associated  with  the  equinoxes,  but  also  exhibit  structures  that  are  unrelated  to 
scintillation  (Figure  [3~3|. 
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Figure  3.7:  Day  of  year  versus  outage  duration  for  Ascension  (a)  and  Ancon  (b) 

in  2001,  and  Parepare  (c)  and  Marak  Parak  (d)  in  2000  with  a  15°  elevation  mask 
applied  and  outages  below  20°  excluded. 
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Table  3.2:  Long  Duration  Outages  at  Marak  Parak,  2000.  The  first  column  of  the 
table  indicates  the  LITC  day  of  year  in  which  the  outage  occurred.  The  duration 
of  the  outage  in  seconds  is  shown  in  the  second  column.  Columns  3-5  indicate  the 
maximum  .Sfi  before,  during  and  after  the  outage  as  reported  by  the  available  GPS 
satellites.  The  final  two  columns  indicate  the  elevation  angle  of  the  affected  satellite 
immediately  before  and  after  the  outage. 

Day  Length  S^ZZ9  Elev  Before  Elev  After 


1 

5820 

0.06 

0.11 

0.04 

34.97 

34.63 

2 

5880 

0.06 

0.08 

0.05 

33.44 

-34.22 

3 

5820 

0.1 

0.08 

0.09 

33.41 

34.2 

4 

5820 

0.14 

0.28 

0.08 

33.39 

34.18 

5 

5820 

0.12 

0.12 

0.08 

33.36 

34.17 

50 

3120 

0.04 

5.83 

0.06 

-40.14 

-18.25 

88 

960 

1.1 

0.34 

9.00 

23.6 

29.27 

89 

660 

0.74 

0.2 

0.94 

20.05 

23.68 

110 

1140 

0.75 

0.71 

0.68 

-21.05 

-18.59 

123 

900 

1.15 

0.17 

0.3 

21.83 

26.0 

285 

2040 

0.2 

0.99 

3.7 

37.54 

45.05 

292 

660 

1.06 

0.31 

0.75 

-20.24 

-18.08 

Table  3.3:  Correlation  Coefficients  for  Max  vs  a2  S 4  showing  the  strong  correlation 
between  these  parameters  before  and  after  an  outage.  The  absence  of  the  satellite 
during  the  outage  degrades  the  correlation,  as  expected. 


S4  Max  and  a2  Outage  Correlation  Coefficients 


Site 

Before  Outage 

During  Outage 

After  Outage 

Ancon 

.94 

.80 

.96 

Ascension 

.95 

.76 

.95 

Marak  Parak 

.95 

.59 

.95 

Parepare 

.94 

.76 

.95 

44 


(a) 


Log  Max  S4 

(c) 


(b) 


(d) 


(e)  (f) 


Figure  3.8:  Relationship  between  maximum  64  and  <r|4  at  Ascension  (left  column) 
and  Ancon  (right  column)  one  minute  before  (a  and  b),  during  (c  and  d)  and  one 
minute  after  (e  and  f)  each  bounded  outage  in  2001.  The  correlation  coefficients  are 
given  in  Table  T3  These  figures  seem  to  suggest  the  geometry  of  GPS  satellites  and 
the  scale  of  EPBs  leave  a  minority  of  satellites  susceptible  to  scintillation,  in  most 
circumstances. 
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Tabic  3.4:  Fraction  of  Satellites  Involved  in  Outages  (a)  &  Percentage  of  Outages 
by  Maximum  S4  (b).  Tabic  3.4  (a)  shows  that  most  outages  involve  a  small  fraction 
of  satellites.  Tabic  3.4  (b)  shows  that  most  simultaneous  outages  occur  during  weak 
scintillation.  These  statisites  exclude  satellites  whose  elevation  is  below  15° 


(a) 


Fraction  of  Available  Satellites  vs.  All  Outages  (%) 


Fraction  of  Satellites 

Ascension 

Ancon 

Parepare 

Marak  Parak 

0 

0 

0.02 

0 

2.42 

0.1 

57.62 

52.86 

46.04 

92.34 

0.2 

23.22 

11.01 

7.69 

4.71 

0.3 

1.85 

0.91 

0.71 

0.03 

0.4 

0.42 

0.82 

0.26 

0.03 

0.5 

0.88 

1.87 

0.39 

0 

0.6 

0.3 

1.35 

0.52 

0.03 

0.7 

0.66 

1.11 

1.53 

0 

0.8 

2.92 

6.31 

8.48 

0 

0.9 

0 

0 

0.39 

0.03 

1 

12.12 

23.73 

33.98 

0.4 

(b) 

Maximum  S4  vs.  Simultaneous  Outages  (%) 


*5*4 

Ascension 

Ancon 

Parepare 

Marak  Parak 

0 

55.69 

55.57 

61.53 

0 

0.1 

29.46 

20.39 

18.1 

50 

0.2 

3.35 

9.92 

5.07 

16.67 

0.3 

2.46 

5.14 

3.38 

16.67 

0.4 

2.12 

3.37 

2.75 

0 

0.5 

1.34 

1.68 

2.36 

0 

0.6 

1.79 

0.94 

1.64 

8.33 

0.7 

0.89 

1.31 

1.5 

0 

0.8 

0.11 

0.65 

0.97 

0 

0.9 

0.33 

0.19 

0.43 

0 

1 

2.46 

0.84 

2.27 

8.33 
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Outages  are  generally  clustered  around  low  elevation  angles,  with  a  second  peak 
between  40-60°  during  ascent  (Figure 

Long  duration  outages  are  primarily  confined  to  elevation  angles  <20°  (Fig¬ 
ure 


3.4). 


After  removing  outages  <20°,  the  majority  of  the  remaining  outages  are  10  min 
or  less  (Figure  [3/7]) . 

Outages  involving  all  satellites  occur  in  as  many  as  ~  34%  of  the  cases,  and  the 
majority  of  these  simultaneous  outages  last  two  minutes  or  less  and  are  associ¬ 
ated  with  5*4  values  below  0.3.  This  suggests  phenomena  other  than  scintillation 


at  work  (Table  pL4[). 

We  said  earlier  in  this  section  ignoring  outages  may  lead  to  under-estimates  in 
the  intensity  and  frequency  of  scintillation.  Consequently,  we  explored  imputing  val¬ 
ues  when  we  were  reasonably  certain  the  data  gap  arose  because  of  scintillation.  The 
observations  above  were  used  to  build  an  algorithm  that  imputed  missing  S4  values 
prior  to  any  further  calculations.  The  decision  to  embark  on  an  imputation  scheme, 
furthermore,  was  made  with  the  understanding  that  any  subsequent  analysis  using 


imputed  data  would  require  modifications  to  standard  statistical  measures  Little  and 


Rubin,  2002  ,  and  that  imputation  has  the  potential  to  distort  the  data.  In  its  final 


form  the  imputation  scheme  worked  as  follows: 


1.  Identify  gaps  in  data  lasting  longer  than  1  min,  but  less  than  10  min.  The  upper 
limit  of  10  min  was  based  on  the  results  after  removing  outages  less  than  20°. 

2.  Outages  lasting  longer  than  10  min  are  ignored. 

3.  If  the  elevation  of  the  missing  satellite  was  last  reported  at  or  below  20°  do  not 
fill  since  this  may  be  a  mask-horizon  problem  or  multipath. 

4.  If  the  elevation  of  the  satellite  before  the  outage  was  greater  than  20  °,  assume 
scintillation  and  impute  a  value  of  8.88  as  a  flag.  Because  we  were  concerned 
only  with  the  number  and  fraction  of  satellites  reporting  £4  above  0.3,  0.5  and 
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0.8,  the  8.88  imputation  value  served  to  place  the  particular  missing  satellite  in 
the  highest  bin,  indicating  severe  scintillation. 

At  the  end  of  our  investigation,  we  attempted  to  confirm  some  solar-cycle, 


seasonal  and  diurnal  behavior  in  the  number  and  fraction  of  satellites  scintillated. 
To  do  this,  we  combined  data  from  various  years,  days  and  times  in  order  to  fill  in 
the  large  scale  gaps.  This  was  a  so  called  “hot-deck”  imputation  scheme  in  which 
missing  values  are  filled  with  values  from  similar  circumstances  (e.g.  same  site,  day 
of  year,  time,  and  solar-cycle  conditions,  but  from  a  different  year).  In  Figure  3.9|  (a) 
each  7-day  week  of  the  year  is  plotted  along  the  abscissa  and  each  site  and  year 
combination  considered  is  plotted  on  the  ordinate.  Black  shading  indicates  at  least  4 
days  of  data  in  a  given  7-day  period  are  available.  White  shading  indicates  areas  where 


data  availability  falls  below  that  threshold.  Figure  3.9  (b)  shows  the  F10.7  cm  solar 
flux  values  and  identifies  the  site/year  combinations  chosen  for  analysis.  A  further 
description  of  the  process  used  to  select  the  data,  and  the  types  of  calculations  carried 


out  is  presented  in  Section  3.3 


3. 3  Data  Processing 

Our  original  objectives  were:  To  determine  the  percent  probability  of  occurrence 
of  S4  exceeding  0.3,  0.5  and  0.8  at  the  50th,  75th  and  95th  percentile  levels  as  a  function 
of  week  and  sunset-relative  time,  and  determine  the  resulting  effect  on  GPS  horizontal 


dilution  of  precision,  HDOP  (see  Subsection  2.1.4).  We  used  Matlab®  to  process  the 
data  as  follows: 


1.  Runs  were  conducted  with  and  without  imputation  for  comparison. 

2.  All  times  were  adjusted  to  be  relative  to  local  sunset  and  “wrapped”  (see  below). 

3.  Seven  days  of  S4  values  were  collected  and  segregated  into  15  minute  “blocks”. 

4.  Counts  were  taken  of  the  number  of  S4  values  exceeding  the  0.3,  0.5  and  0.8 
thresholds  for  a  given  time  block  for  each  of  the  seven  days. 
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Weekly  Data  Availability  (>  4  Days  Available) 


c 

o 

H — ' 

CO 

o 

o 

_J 

T3 

C 

CO 

■> _ 

CO 

05 
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2004  parepare 
2004  marakparak 
2004  ascension 
2004  ancon 
2003  parepare 
2003  marakparak 
2003  ascension 
2003  ancon 
2002  parepare 
2002  marakparak 
2002  ascension 
2002  ancon 
2001  parepare 
2001  marakparak 
2001  ascension 
2001  ancon 
2000  parepare 
2000  marakparak 
2000  ascension 
2000  ancon 
1 999  parepare 
1999  marakparak 
1999  ascension 
1 999  ancon 
1 998  parepare 
1 998  marakparak 
1998  ascension 
1 998  ancon 
1 997  parepare 
1 997  marakparak 
1997  ascension 
1 997  ancon 


10  13  16  19  22  25  28  31  34  37  40  43  46  49  52 
Week  of  Year 


(a) 


ISES  Solar  Cycle  FI 0.7cm  Radio  Flux  Progression 

Data  Through  31  Dec  06 


- Smoothed  Monthly  Volues  — * —  Monthly  Values 

Updated  2007  Jan  3  NOAA/SEC  Boulder, CO  USA 


(b) 


Figure  3.9:  Large  scale  data  gaps  (a)  and  solar  cycle  relationships  (b).  White 

regions  in  (a)  indicate  missing  data,  black  indicates  at  least  4  days  of  data  are 
available  in  the  corresponding  week.  In  (b),  the  years  chosen  to  be  grouped  together 
for  each  site  are  shown  in  relation  to  the  solar  activity  as  indicated  by  the  F10.7  cm 
solar  flux. 
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5.  For  each  £4  threshold,  the  data  was  sorted  and  the  50th,  75th  and  95th  percentiles 
calculated. 

6.  Resulting  values  were  plotted  by  week  and  year. 


Since  we  were  concerned  with  nocturnal  scintillation  (section  2.2.2),  the  times  were 
adjusted  from  UTC  to  local  sunset  for  easy  comparison  between  sites.  For  some  loca¬ 
tions,  data  spanned  UTC  midnight;  in  other  words,  a  single  night’s  scintillation  was 
spread  over  two  UTC  days.  Consequently  there  was  a  portion  of  data  in  the  early 
morning  hours  on  a  given  UTC  day,  a  large  gap  during  daylight,  then  more  data  be¬ 
ginning  in  the  UTC  evening.  For  computational  efficiency,  rather  than  concatenating 
data  records,  the  data  prior  to  sunrise  on  a  given  UTC  day  was  “wrapped”  around  to 
the  end  of  that  UTC  day  by  adding  86400  s  to  the  sunset-relative  time  stamps.  For 
the  computations  above,  this  adjustment  was  transparent. 


Collecting  seven  days  of  data  in  15  minute  “blocks”  provided  a  maximum  of 
105  observations  per  satellite  per  block,  assuming  all  observations  were  recorded  or 
imputed.  We  also  required  420  observations  per  block,  arriving  at  this  threshold  by 
assuming  at  least  four  satellites  were  reporting  over  the  seven  day  period;  4  x  105  = 
420.  When  the  total  number  of  observations  fell  below  420,  calculations  were  not 
performed  and  the  block  was  skipped. 


HDOP  was  calculated  at  one  minute  intervals  for  the  same  data.  Initially, 
HDOP  was  calculated  using  all  available  satellites,  regardless  of  £4  value.  Then  satel¬ 
lites  with  S4  values  exceeding  the  thresholds  above  were  incrementally  removed  and 
new  HDOPs  calculated.  Recall  from  Chapter  [TT]  in  order  to  calculate  HDOP,  the 
satellite  positions  are  required.  However,  for  the  imputed  data,  no  position  informa¬ 
tion  was  available.  Consequently,  almanac  datc^j  were  used  in  the  following  manner 
to  calculate  HDOP: 


1.  Almanac  data  were  compiled  for  the  particular  day  and  times  for  which  obser¬ 
vational  data  was  available,  using  code  provided  by  Dr.  John  Raquet. 


4subsection 


2.1.2 


50 


2.  The  H  matrix  ,Equation|2.7[  was  loaded  with  the  positions  of  all  satellites  above 
a  20°  mask. 


3.  If  almanac  output  indicated  a  visible  satellite  that  did  not  appear  in  the  data, 
that  satellite  was  assigned  an  S4  value  of  999.99. 

4.  HDOP  was  calculated  for  the  following  thresholds:  All  satellites  regardless  of 
S4  ,  S4  <  0.8,  0.5  and  0.3.  If  less  than  four  satellites  were  available,  an  HDOP 
value  of  -1  was  assigned,  indicating  the  inability  to  calculate  HDOP. 

5.  All  times  were  adjusted  to  be  sunset  relative  and  data  was  “wrapped  around”  as 
before.  The  effect  of  this  “wrapping”  will  be  apparent  when  one  minute  HDOP 
values  are  plotted. 


At  the  end  of  our  investigation,  we  also  created  summary  plots  of  diurnal,  sea¬ 


sonal  and  solar-cycle  behavior  using  hot-deck  imputation  described  in  Section  3.2 


The  same  percentiles  and  thresholds  were  used,  but  no  one-minute  imputation  was 
carried  out. 


These  procedures  were  scripted  in  Matlab®  and  applied  to  the  data.  Both  the 
creation  of  the  statistics  and  figure  generation  were  automated  to  allow  an  investigator 
to  choose  any  combination  of  sites  and  years  to  analyze.  In  the  next  chapter,  we 
present  the  results  of  the  investigation  and  corresponding  images. 
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Table  3.5:  Summary  Plot  Data  Combinations  showing  years  used  and  grouping  of 
one-minute  data.  For  the  solar  cycle  comparison,  individual  years  were  analyzed.  For 
the  other  comparisons,  the  data  from  all  the  years  shown  was  combined.  One-minute 
data  was  grouped  into  15min  blocks  or  simply  combined.  For  the  seasonal  comparison, 
all  weeks  were  treated  individually.  For  the  solar  cycle  and  diurnal  comparison,  weeks 
11-13  and  37-39  (centered  on  the  equinoxes)  were  used. 


Solar  Cycle 

Seasonal 

Diurnal 

Ascension 

2000 

1999 

1999 

2004 

2002 

2003 

2002 

Ancon 

1998 

1999 

1999 

2001 

2002 

2003 

2002 

Parepare 

1998 

1998 

1999 

2000 

1999 

2002 

2002 

Marak  Parak 

1999 

1998 

1999 

2001 

1999 

2002 

2002 

One  Minute  Data 

15-min  block 

all 

15-min  block 

Weeks 

11-13  &  37-39 

Individual 

11-13  &  37-39 
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IV.  Results  and  Analysis 


Now  we  examine  the  statistics  generated  using  the  procedures  presented  in  Chap¬ 
ter  El  We  investigate  the  solar  cycle,  seasonal  and  diurnal  behavior  described 
in  Chapter  [TTJ  This  analysis  is  carried  out  for  both  raw  and  imputed  data.  One- 
minute  HDOP  results  are  evaluated.  Finally,  we  scrutinize  the  imputation  results. 
Before  embarking  on  this  journey,  however,  we  briefly  present  some  examples  of  our 
work  and  provide  a  commentary  on  interpreting  the  output. 


In  Section  3.3  we  said  we  were  going  to  examine  the  percent  probability  of  oc¬ 
currence  of  exceeding  three  S±  thresholds  at  three  percentile  levels.  The  S4  thresholds 
chosen  were  0.3,  0.5,  and  0.8  and  represent  weak-moderate,  moderate-strong  and  ex¬ 
treme  levels  of  scintillation.  The  50i/!  (also  called  the  median),  75th  (also  called  the 
upper  quartile),  and  95th  percentiles  were  chosen.  Wilks  |1995|  provides  the  following 
interpretation  of  percentiles: 


A  sample  quantile  [percentile]  qp  is  a  number  having  the  same  units  as  the 
data,  which  exceeds  that  proportion  of  the  data  given  by  the  subscript  p 
with  0  <  p  <  1.  The  sample  quantile  qp  can  be  interpreted  approximately 
as  the  data  value  exceeding  a  randomly  chosen  member  of  the  data  set, 
with  probability  p.  Equivalently,,  the  sample  quantile  qp  would  be  regarded 
as  the  [p  x  100]t/l  percentile  of  the  data  set. 


Stated  another  way,  there  is  a  (1  —  p)  x  100  percent  chance  that  a  randomly 
sampled  value  would  exceed  that  found  at  the  [p  x  100]t/l  percentile  of  the  data  set. 


Figure  [47T|  provides  an  example  for  discussion.  Figure  [47T|  is  a  plot  of  the  percentage 
of  available  satellites  (indicated  on  the  ordinate)  being  scintillated  with  S4  >0.5  as 
a  function  of  hours  after  local  sunset  (indicated  on  the  abscissa)  for  Ascension  in 
2001  during  week  11.  The  data  has  been  sorted  into  15-minute  bins  as  described 
in  Section  |3.3|  and  no  imputation  has  been  carried  out.  Three  lines  represent  the 


percentage  of  satellites  at  the  50^  (heavy  solid  line),  75th  (dash-dot  line)  and  95th 
(thin  solid  line)  percentiles.  Four  hours  after  local  sunset  we  can  see  that  at  the 
50^  percentile,  slightly  under  35%  of  the  available  satellites  can  be  expected  to  have 
an  S4  >0.5.  At  the  75th  percentile,  about  4.5  h  after  local  sunset,  the  number  of 
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satellites  reporting  S4  >0.5  is  60%  while  at  the  95™  percentile  the  value  climbs  to 


approximately  89%.  This  is  expected  since,  as  noted  in  Section  322  high  values  of 
S4  typically  affect  a  small  percentage  of  satellites. 

In  addition  to  weekly  plots,  we  also  constructed  yearly  summary  plots.  Fig- 
4.3  is  an  example.  In  Figure  [4~3|  the  abscissa  indicates  the  sunset  relative  time  as 


ure 


before,  however  the  ordinate  now  shows  week-number.  The  colors  represent  the  vari¬ 
ous  fractions  of  satellites  affected  by  scintillation  at  the  particular  S4  and  percentile 
combination.  In  this  figure,  the  median  (a),  upper  quartile  (b)  and  95th  percentile  (c) 
for  S4  >0.3  are  displayed.  The  white  areas  indicate  insufficient  data.  As  before,  no 
imputation  was  carried  out.  The  striations  that  run  from  the  upper  left  to  lower  right 
of  the  figure  are  a  result  of  the  precession  of  the  GPS  satellites.  Each  day  the  satellite 
arrives  four  minutes  earlier  than  the  day  before;  the  number  of  available  satellites 
change  in  regular  patterns  throughout  the  year.  These  figures  show  the  same  increase 
in  the  fraction  of  affected  satellites  with  increasing  percentile  level.  I11  short,  the 
combination  of  high  percentile  and  low  £4  threshold  will  produce  the  most  dramatic 
results.  Consequently,  the  bulk  of  subsequent  figures  will  have  an  S4  threshold  of  0.3 
at  the  95th  percentile. 


4-1  Solar  Cycle  Comparison 


Recall  from  Table  2.2.2  the  production  of  EPBs  and  resulting  scintillation  tends 


to  mirror  the  solar  cycle.  In  Figure  4.4  shows  the  solar  cycle  sensitivity  for  As¬ 
cension  (a),  Ancon  (b),  Parepare  (c)  and  Marak  Parak  (d).  These  plots  depict  the 
fraction  of  available  satellites  experiencing  S4  >0.3  at  the  95th  percentile  as  a  func¬ 
tion  of  hours  after  local  sunset.  For  each  site,  years  for  comparison  were  selected  near 


solar  maximum  and  solar  minimum  as  indicated  in  Figure  3.9  (b).  For  all  the  sites  in 


Figure  [4~4l  the  influence  of  the  solar  cycle  is  apparent.  The  effect  is  most  pronounced 
at  the  Appleton  Anomaly  crest  stations  of  Ascension  (a)  and  Parepare  (c),  although 
the  fraction  of  affected  satellites  jumps  as  much  as  30%  at  the  geomagnetic  equator 
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Figure  4.1:  Ascension  2001,  Week  11,  Percentage  of  Satellites  with  S4  >  0.5  at  the 
50ih  (heavy  solid  line),  75th  (dash-dot  line)  and  95th  (thin  solid  line)  percentiles.  This 
figure  was  created  using  raw  data. 


stations.  This  is  consistent  with  the  behavior  presented  in  Table  2.2.2  which  noted 
the  bubble  activity  was  highest  under  the  anomaly  crests. 

We  also  examined  the  impact  on  the  results  using  our  imputation  scheme  de¬ 
scribed  in  Section [T2  Figure  [475  shows  the  comparison  for  Parepare  during  1998  (a) 
&  (c)  and  2000  (b)  &  (d).  In  this  example,  (a)  and  (b)  represent  the  raw  data, 
while  data  gaps  in  (c)  and  (d)  were  imputed.  The  abscissa  indicates  the  hours  after 
local  sunset;  the  ordinate  indicates  the  week  number.  The  color  scale  represents  the 
fraction  of  available  satellites  with  S4  >0.3  at  the  95th  percentile.  Again  the  solar 
cycle  influence  is  apparent,  with  much  higher  values  appearing  in  2000  during  the  so¬ 
lar  maximum,  particularly  during  the  spring  equinox.  The  impact  of  the  imputation 
scheme  is  less  obvious,  although  close  inspection  will  reveal  enhanced  values  between 


weeks  30  and  35  and  during  week  40  for  the  1998  data  (Figure  4.5  (a)  and  (c)).  The 


impact  of  imputation  will  be  discussed  in  greater  detail  in  Section  4.5 


55 


(a) 


(b) 


Figure  4.2:  Ascension  2001,  Week  11,  Percentage  of  Satellites  with  S4  >  0.3  (a) 
and  S4  >  0.8  (b).  These  figures  were  created  using  raw  data. 


1  2  3  4  5  6  7 

Hours  After  Local  Sunset 
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Figure  4.3:  2001  Summary  for  Ascension,  Ratio  of  satellites  with  S4  >0.3  to  all 

available,  50ih  (a),  75th  (b),  and  95th  (c)  Percentile.  The  abscissa  indicates  the  hours 
after  local  sunset;  the  ordinate  indicates  the  week  number.  Striations  running  from 
upper  left  to  lower  right  are  an  artifact  of  GPS  precession:  the  number  of  satellites 
available  changes  in  regular  patterns  throughout  the  year.  These  figures  were  created 
using  raw  data. 
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Figure  4.4:  Solar  cycle  sensitivity  for  Ascension  (a),  Ancon  (b),  Parepare  (c)  and 

Marak  Parak  (d).  These  plots  depict  the  fraction  of  available  satellites  with  S4  >0.3 
at  the  95th  percentile  as  a  function  of  hours  after  local  sunset.  The  abscissa  indicates 
the  hours  after  local  sunset;  the  ordinate  indicates  the  week  number.  These  figures 
were  created  using  raw  data 


Figure  4.5:  Solar  Cycle  Comparison  Parepare  for  1998  (a)  &  (c)  and  2000  (b)  &  (d). 
Raw  data  was  used  for  (a)  and  (b)  while  (c)  and  (d)  were  filled  using  the  imputation 
scheme  described  in  Section  13.21 
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4-  2  Seasonal  Comparison 


Figure  4.6  shows  the  results  of  a  seasonal  comparison  using  data  from  the  years 


identified  in  Table  3J3  Like  the  earlier  figures,  the  values  on  the  ordinate  represent  the 
fraction  of  available  satellites  experiencing  £4  levels  above  0.3  at  the  95th  percentile 
at  the  four  sites.  The  abscissa,  however,  now  indicates  the  week  number.  Both  the 
data  (markers)  and  the  smoothed  results  (solid  lines)  were  plotted.  The  smoothed 
results  were  created  by  applying  a  five-week  boxcar-average  to  the  raw  results.  Other 
than  the  “hot  deck”  imputation  carried  out  by  combining  the  different  years,  no  other 
imputation  was  performed. 


In  Figure  |4.6[  the  equinox  activity  peaks  described  in  Table  |2.2.2|  are  evident, 
although  the  onset  and  cessation  of  the  enhanced  scintillation  is  not  precisely  centered 
on  the  equinox.  Ascension’s  activity  seemed  to  reach  its  zenith  on  the  equinox,  while 
Parepare’s  was  delayed  by  2-3  weeks.  At  the  geomagnetic  equator,  Ancon’s  maxima 
appeared  to  occur  around  week  6,  while  Marak  Parak  did  not  reach  its  spring  peak 
until  week  21.  The  response  was  more  uniform  around  the  autumnal  equinox  with 
all  sites  showing  a  peak  in  activity.  Unlike  the  spring,  Marak  Parak  showed  a  well- 
defined  maximum  and  was  the  closest  to  the  autumnal  equinox.  The  maximum  levels 
at  Ascension  and  Parepare  were  almost  simultaneous,  while  Ancon  lagged  behind. 
The  amplitude  of  the  seasonal  cycle  was  higher  for  the  anomaly  crest  stations,  but  all 
reflected  the  solstitial  lull. 


Figure  4/7  is  the  seasonal  comparison  for  Ascension  showing  the  fraction  of 
satellites  with  S4  >0.3  at  the  95th  percentile.  The  abscissa  indicates  the  hours  after 
local  sunset;  the  ordinate  indicates  the  week  number.  Figure |4~7|  (a)  shows  the  results 
for  the  raw  version,  while  (b)  shows  the  results  for  the  imputed  version.  The  equinox 
peaks  mentioned  in  the  preceding  paragraph  are  visible  in  both  the  imputed  and  raw 
versions,  although  the  fall  equinox  values  are  somewhat  enhanced  in  the  imputed 
version.  Here  S4  >0.3  affects  all  available  satellites  and  lasts  for  over  an  hour  during 
weeks  14,  42-45. 
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Week 


Figure  4.6:  Seasonal  Comparison  for  Ascension,  Ancon,  Parepare,  and  Marak 

Parak;  fraction  of  satellites  with  S4  >0.3  at  the  95th  percentile.  Solid  lines  repre¬ 
sent  the  reuslts  when  a  5-week  filter  was  applied  to  the  results.  The  original  data 
points  are  indicated  by  the  markers. 
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Figure  4.7:  Seasonal  Comparison  for  Ascension  showing  the  fraction  of  satellites 

with  S4  >0.3  at  the  95th  percentile.  The  abscissa  indicates  the  hours  after  local 
sunset;  the  ordinate  indicates  the  week  number.  Figure  4.7  (a)  shows  the  results 
using  raw  data,  while  (b)  shows  the  results  using  the  imputation  scheme  described  in 
Section  13.21 
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4-3  Diurnal  Comparison 


Recall  in  Table  |2.2.2|  the  diurnal  cycle  was  characterized  by  the  onset  of  EPB 
and  scintillation  activity  about  an  hour  after  local  sunset,  a  peak  around  2200L  and 


dissipation  by  sunrise.  Figure  4.8  shows  the  diurnal  behavior  for  Ascension  and  An¬ 


con  (a)  and  Parepare  and  Marak  Parak  (b),  using  data  from  the  years  identified  in 
Table  [375]  These  ordinate  indicates  the  fraction  of  available  satellites  with  S4  >0.3  at 
the  95th  percentile  as  a  function  of  hours  after  local  sunset  indicated  on  the  abscissa. 


de  La  Beaujardiere  et  al.  2004  notes  that  begins  at  the  anomaly  crest  stations  sixty 


to  ninety  minutes  later  than  the  equatorial  stations.  She  attributes  this  to  . .  a  finite 
upwclling  speed  of  the  irregularities  at  the  magnetic  equator.”  This  behavior  is  most 


evident  in  the  Atlantic  &  Americas  longitudinal  sector  stations,  Figure  4.8  (a),  where 


the  ramp-up  of  activity  occurs  about  three  hours  after  local  sunset  at  Ascension, 
while  the  increase  at  Ancon  occurs  at  the  two-hour  point.  In  the  Pacific  longitudinal 
sector,  however,  the  lag  appears  to  be  15  to  30  minutes.  This  may  be  attributed  to 
the  relative  differences  in  location.  Ascension  and  Ancon  are  widely  separated  by 
latitude  and  longitude,  while  Parepare  and  Marak  Parak  are  at  essentially  the  same 


longitude  (see  Figure  3.1  Also,  Table  3.1  shows  that  Ancon  is  closer  to  the  magnetic 


dip  equator  than  Marak  Parak. 


Figure  T9  shows  the  diurnal  cycle  at  Parepare  for  week  13  during  2000  for 
the  median,  upper  quartile  and  95th  percentile.  The  ordinate  is  now  the  number  of 
satellites  with  S4  >0.3,  and  the  abscissa,  as  before,  is  the  hours  from  local  sunset. 


Figure  4.9  (a)  shows  the  raw  results,  while  (b)  shows  the  imputed  results.  The  diurnal 
cycle  is  again  well  represented,  and  the  influence  of  imputation  is  more  apparent.  At 
the  median  level,  there  are  times  when  the  number  of  satellites  affected  by  S4  >0.3 
rises  to  three  in  the  imputed  case,  from  two  in  the  raw  version. 


4-4  Annual  Horizontal  Dilution  of  Precision  (HD OP)  Summaries 

HDOP  calculations  were  made  using  the  technique  described  in  Section 


Figure  |4~T0  shows  the  results  for  Ascension  in  2001  using  raw  one-minute  data.  Fig- 
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Figure  4.8:  The  diurnal  behavior  for  Ascension  and  Ancon  (a)  and  Parepare  and 
Marak  Parak  (b),  using  data  from  the  years  identified  in  Table  3.5  These  ordinate 
indicates  the  fraction  of  available  satellites  with  S4  >0.3  at  the  95i/!  percentile  as  a 
function  of  hours  after  local  sunset,  which  is  indicated  on  the  abscissa. 
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Figure  4.9:  Diurnal  Comparison,  Parepare  week  13,  2000,  raw  (a)  and  imputed  (b) 
results.  The  diurnal  cycle  at  Parepare  for  the  median,  upper  quartilc  and  95th  per¬ 
centile.  The  ordinate  is  now  the  number  of  satellites  experiencing  S±  >0.3,  and  the 
abscissa,  as  before,  is  the  hours  from  local  sunset. 
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ure  4.10  (a)  shows  HDOP  obtained  assuming  no  scintillation  but  using  only  satellites 
shown  in  the  data,  as  opposed  to  all  satellites  indicated  by  the  almanac.  Consequently, 


some  periods  where  HDOP  could  not  be  determined  are  still  evident.  Figures  4.10  (b) 
through  (d)  show  the  results  after  removing  all  satellites  with  an  S 4  >0.3,  0.5  and 
0.8  respectively.  The  ordinate  in  this  case  is  the  day  of  year,  and  the  abscissa  is  the 
time  after  local  sunset.  The  color  bar  indicates  the  HDOP  value.  Recall  from  Chap¬ 
ter  [TT|  that  HDOP  =  1  is  considered  perfect.  The  white  spaces  indicate  missing  data 
while  the  gray  shading  indicates  the  inability  to  calculate  HDOP  because  less  than 
four  satellites  were  available  below  the  particular  S4  threshold.  As  the  S4  threshold 
increases  from  (c)  to  (d)  and  the  number  of  available  satellites  increases,  notice  how 
the  HDOP  values  drop  and  amount  of  gray  shading  diminishes.  In  Figure  4.11  we 


have  plotted  the  relative  difference  in  HDOP  between  the  S4  >0.3  threshold  case 
Figure  4.10  (b)  and  the  case  in  which  no  scintillation  was  occurring  Figure  4.10|  (a). 
The  largest  relative  changes  (i.e.  1)  are  coincident  with  the  areas  of  gray  shading  in 


Figure  4.10  (b)  where  no  solution  was  available.  Less  dramatic  yet  still  substantial 


changes  are  peppered  around  these  areas  in  Figure  4.11 


4-5  Imputation  Performance 


Figures  4.12  (a)  and  (b)  depict  the  occurrence  of  differences  in  the  fraction  of 


scintillated  satellites  between  imputed  and  raw  results  for  Ascension  in  2001.  The 


abscissa  is  hours  past  sunset  and  the  ordinate  is  the  week.  Figure  4.12  (a)  shows  the 


occurrence  of  differences  between  the  imputed  and  raw  results  at  the  median  for  the 
fraction  of  satellites  with  £4  >0.8,  while  (b)  shows  the  95th  percentile  for  the  fraction 
of  satellites  with  £4  >0.3.  The  changes  are  confined  to  the  equinox  peaks  in  (a), 
while  in  (b)  the  seasonal  clusters  are  ill-defined  and  changes  begin  to  appear  in  the 
off-season  as  well.  Similar  results  were  obtained  for  Parepare  (not  shown),  although 
the  correlation  are  not  as  compact.  These  comparisons  were  made  during  the  solar 
maximum.  During  solar  minimum  (not  shown),  the  patterns  are  more  diffuse,  with 
fewer  difference  between  raw  and  imputed  percentages.  Ascension  and  Parepare  are 
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Figure  4.10:  HDOP  Comparison:  HDOP  results  without  any  satellites  removed  for 
scintillation  is  shown  in  Figure  4.10  (a)  while  Figure  4.10  (b)  through  (d)  show  the 
results  after  removing  all  satellites  with  an  S4  >0.3,  0.5  and  0.8  respectively.  The 
ordinate  in  this  case  is  the  day  of  year,  and  the  abscissa  is  the  time  after  local  sunset. 
The  color  bar  indicates  the  HDOP  value  (recall  from  Chapter  [TT]  that  a  value  of  1 
is  considered  perfect.)  Gray  areas  indicate  times  when  less  than  four  satellites  were 
available,  and  hence  no  navigation  solution,  or  HDOP  calculation,  was  available. 
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4.11 


is  the  relative  difference  in  HDOP  between  the  Sh  >0.3  threshold 


Figure  4.10  (b)  and  the  case  in  which  no  scintillation  was  occurring  Figure  4.10  (a). 


both  situated  under  the  anomaly  crest.  We  examined  Ancon  and  Marak  Parak  (not 
shown),  closer  to  the  geomagnetic  equator.  Here  the  seasonal  correlation  was  much 
harder  to  discern,  particularly  during  solar  minimum.  Based  on  an  examination  of 


plots  like  those  shown  in  Figure  |4.12[  the  performance  of  the  imputation  scheme 
appears  to  deteriorate  as  percentile  increases  and  S4  and  solar  activity  decrease.  To 
gain  a  better  understanding  of  the  magnitude  of  the  difference  between  the  imputed 
and  raw  results,  we  performed  a  simple  hypothesis  test.  For  each  of  the  32  15- 
min  blocks,  for  each  week  at  a  given  location,  year,  percentile  and  S4  threshold, 
we  calculated  the  difference  di  between  the  imputed  and  raw  results.  If  there  was 
110  significant  difference  between  the  raw  and  imputed  results,  the  mean  of  their 
differences,  should  be  zero.  If  Ha  7^  0,  the  imputation  scheme  had  a  statistically 
significant  impact.  This  can  be  expressed  as  a  null  and  alternative  hypothesis  as: 


Ho  •  l^d  0 
Ha-Hd.  7^  0 


We  chose  a  Z-test  on  the  mean  of  the  differences  Hd  at  the  a  =  0.05  significance 
level.  In  order  to  apply  the  test,  we  had  to  ensure  the  samples  were  independent.  We 
examined  the  lag-1  autocorrelation  for  the  differences  and  found  the  overwhelming 
majority  (>  84%)  of  the  samples  were  independent.  Reasonably  reassured  by  the 
autocorrelation  results,  we  applied  the  Z-test  to  all  sites,  percentiles  and  S4  levels 
during  solar  minimum  and  solar  maximum.  In  all  cases  we  were  able  to  reject  the 
null  hypothesis  at  the  a  =  0.05  significance  level,  meaning  there  was  a  statistically 
significant  difference  between  the  imputed  and  raw  versions.  The  results  are  tabulated 
in  Appendix  [B]  and  show  the  maximum  of  the  95%  confidence  interval  for  the  mean 
reached  14%,  however  values  were  typically  less  than  10%.  A  scatter  plot  of  the 
hypothesis  test  results  showed  the  highest  mean  differences  were  typically  near  the 
equinoxes,  as  expected.  Our  analysis  suggests  that  the  imputation  scheme  can  produce 
statistically  significant,  albeit  small,  differences  in  the  fraction  of  satellites  scintillated. 


In  order  to  strengthen  the  imputation  scheme  and  remove  anomalous  off-season 
results,  we  added  a  lower  limit  of  0.3  for  minimum  S4  in  the  minute  before  an  out¬ 
age  for  a  particular  PRN  and  relaxed  the  duration  limit  to  20  min  from  the  10  min 


limit  originally  described  in  Section  |3.2|  We  added  the  S4  threshold  to  ensure  the 
satellite  was  at  least  experiencing  weak  scintillation  before  an  outage.  We  expanded 
the  time  limit  to  20  min  because  a  substantial  number  of  the  outages  of  this  duration 


were  present  even  after  removing  the  outages  below  20°  (see  e.g.  Figure  3.7  (a)). 
This  refinement  produced  the  results  shown  in  Figure  |4.13|  (a)  and  (b).  Note  the 


disappearance  of  off-season  (solstitial)  changes  in  Figure  4.13  (b)  when  compared  to 
Figure  4.12|  (b).  The  number  of  time-blocks  that  increased  from  the  raw  version  also 
rose  as  seen  in  the  appearance  of  additional  +  symbols  around  the  equinox  peaks 
in  both  Figure  4.13|  (a)  and  (b).  We  conducted  a  Z-test  on  the  means  of  the  dif¬ 
ferences  as  we  did  for  the  original  scheme.  The  confidence  interval  limits  for  the 
mean  of  the  differences  were  generally  lower  for  this  imputation  scheme  than  for  the 
original  for  the  solar  minimum  years  and  the  geomagnetic  equator  locations  (Ancon 
and  Marak  Parak).  The  confidence  interval  limits  generally  increased  during  solar 
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Figure  4.12:  Difference  between  Imputed  and  Raw  data  at  Ascension  in  2001  using 

The  ’+’  indicates  an  increase 


imputation  criteria  originally  described  in  Section  |3.2 
over  the  raw  data 


Figure  4.12 


i  shows  the  difference  at  the  median  for  S4  >0.8 
while  Figure  4.12  (b)  shows  the  result  at  the  95th  percentile  for  while  for  S4  >0.3. 
As  in  earlier  figures,  the  abscissa  indicates  the  sunset-relative  time  while  the  ordinate 
indicates  the  week  of  the  year. 
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maximum  at  the  anomaly  crest  locations  (Ascension  and  Parepare).  Additionally, 
the  confidence  interval  peaks  were  more  concentrated  around  the  equinoxes  than  in 
the  original  scheme.  Appendix  [B]  contains  the  comparison.  These  results  suggest  in¬ 
cluding  an  5*4  threshold  and  loosening  the  duration  constraint  produces  more  realistic 
output,  opening  the  door  for  further  refinement  of  the  imputation  scheme.  Finally, 
we  increased  the  S4  threshold  to  0.8  for  a  particular  satellite  prior  to  an  outage,  guar¬ 
anteeing  that  any  satellite  for  which  missing  values  were  imputed  was  experiencing 
severe-extreme  scintillation  in  the  minute  prior  to  being  lost.  The  results  (not  shown) 
were  very  similar  to  those  using  the  0.3  thresholds,  although  the  Z-test  confidence 
intervals  were  in  some  cases  1-2%  lower.  This  further  increased  our  confidence  in  the 
modified  imputation  scheme. 


4-6  Summary 


Thus  we  have  been  able  to  continue  the  work  begun  by  Thomas  et  al.  2001  ,  ex¬ 


panding  it  to  include  stations  in  the  Atlantic  and  Americas  sector,  as  well  as  exploring 
solar  cycle  influences.  We  were  able  to  observe  the  behavior  described  in  Chapter  [hi 
in  most  cases  with  the  following  exceptions: 


•  The  seasonal  behavior  during  the  spring  at  Marak  Parak  was  peculiar,  with  the 
peak  activity  delayed  until  approximately  week  21. 

•  In  general  the  seasonal  peaks  were  near,  but  not  necessarily  centered  on  the 
equinoxes  at  the  four  sites  studied. 

•  The  diurnal  onset  of  scintillation  conditions  at  the  anomaly  crest  site  at  Parepare 
appeared  to  lag  Marak  Parak  by  only  15-30  min,  rather  than  the  60-90  min  we 
expected.  The  delay  was  about  an  hour  in  the  Atlantic  sector.  The  variation  is 
most  likely  caused  by  the  relative  differences  in  location. 

The  HDOP  plots  clearly  show  the  impact  of  scintillation  on  GPS  precision.  As 
satellites  are  removed  based  on  S4  thresholds,  the  HDOP  increases  and  in  extreme 
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Figure  4.13:  Difference  between  Imputed  and  Raw  data  at  Ascension  in  2001  after 
adding  a  lower  limit  of  0.3  for  maximum  S4  in  the  minute  before  an  outage  and 
relaxing  the  duration  criteria  to  20  min  from  the  10  min  criteria  originally  described 


in  Section  3.2  The  ’+’  indicates  an  increase  over  the  raw  data.  Figure  4.13  (a)  shows 


the  difference  at  the  median  for  S4  >0.8  while  Figure  4.13  (b)  shows  the  result  at  the 
95 th  percentile  for  while  for  S4  >0.3.  As  in  earlier  figures,  the  abscissa  indicates  the 
sunset-relative  time  while  the  ordinate  indicates  the  week  of  the  year. 
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cases,  a  navigation  solution  becomes  impossible  because  less  than  four  satellites  re¬ 
main.  This  result  has  serious  ramifications  for  GPS  useres. 

Finally,  we  saw  that  the  imputation  scheme  increases  the  observed  fraction  of 
scintillated  satellites  at  a  typically  small  (<  10%)  yet  statistically  significant  level  and 
that  the  scheme  can  be  improved  by  including  an  S4  threshold  and  relaxing  the  dura¬ 
tion  constraint.  All  that  remains  for  this  work  is  a  summation  and  recommendations 
for  the  future,  which  are  presented  in  the  next  chapter. 
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V.  Summary  and  Conclusions 


WE  have  extended  the  work  of  Thomas  et  al.  1 2001 1  incorporating  new  sites 
and  years.  Additionally,  a  simple  ad-hoc  imputation  technique  was  devel¬ 
oped,  tested  and  refined.  The  synoptic  behavior  that  emerged  from  our  analysis  was 
consistent  with  the  findings  of  other  researchers.  Choosing  to  study  S4  values  col¬ 
lected  during  the  solar  maximum  and  at  locations  beneath  the  Appleton  Anomaly 
crest  resulted  in  some  eye-opening  images  and  statistics.  That  GPS  could  be  unavail¬ 
able  to  a  military  user  for  an  hour  or  more  is  significant,  particularly  if  operation 
is  planned  without  accounting  for  scintillation.  Even  if  GPS  is  available,  precision 
may  be  severely  degraded.  It  is  important  to  remember,  also,  that  the  data  for  this 
research  was  collected  with  fixed  receivers  designed  to  track  through  moderate  severe 
scintillation.  Combat  seldom  provides  such  a  favorable  environment  for  either  the 


equipment  or  the  GPS  user  |  Thomas  et  al.  2001 


We  had  the  advantage  of  comparing  scintillation  activity  across  the  solar  cycle. 
As  anticipated,  there  was  much  more  activity  during  the  solar  maximum  at  all  sta¬ 
tions,  but  again  particularly  at  the  anomaly  crests  in  both  longitude  sectors.  During 
the  years  closer  to  the  solar  minimum,  activity  was  diminished,  but  still  significant. 
We  were  also  able  to  observe  some  departures  from  anticipated  behavior  at  seasonal 
and  diurnal  time  scales.  We  also  discovered  in  the  diurnal  comparisons,  a  longitudinal 
difference  in  anomaly  crest  onset  delay  when  compared  to  equatorial  stations.  The 
HDOP  information  we  produced  provides  a  bridge  from  relatively  esoteric  S4  data  to 
operationally  applicable  dilution  of  precision  information.  This  is  perhaps  the  most 
valuable  portion  of  our  research  for  the  warfighter. 


In  an  effort  to  gain  a  better  estimate  of  scintillation  and  its  operational  im¬ 
pacts,  we  attempted  an  imputation  scheme.  While  enabling  us  to  “fill  in  the  holes” 
and  attempt  to  capture  the  true  magnitude  of  the  scintillation  problem,  the  algorithm 
appeared  over  zealous  when  confronting  low  £4  thresholds  and  high  percentile  levels. 
In  light  of  these  findings,  we  included  a  minimum  £4  threshold  to  further  separate 
scintillation  induced  outages  from  non-scintillation  outages.  We  also  relaxed  the  out- 
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age  duration  limit  to  20  min  to  capture  more  outages.  The  results  were  promising, 
producing  patterns  of  change  over  raw  data  that  were  more  realistic.  Even  after  ap¬ 
plying  a  very  stringent  minimum  S4  threshold  of  0.8,  the  outcome  was  consistent, 
increasing  onr  confidence  in  the  scheme. 

Finally,  to  our  knowledge  this  was  the  Erst  exploration  of  the  patterns  of  “miss¬ 
ingness”  occurring  in  the  64  data.  We  were  able  to  turn  missing  data  into  “miss¬ 
ingness”  data  by  examining  the  circumstances  surrounding  the  gaps.  We  believe  the 
outages  may  represent  another  way  of  investigating  scintillation,  a  means  toward  im¬ 
proving  imputation  schemes,  and  a  potential  climatological  aid  for  planners  to  apply 
operationally.  This  naturally  leads  us  to  a  discussion  of  future  work. 

5. 1  Future  Work 

If  further  statistics  are  to  be  generated  using  imputation,  more  robust  methods 
of  dealing  with  missing  data  are  needed.  We  introduced  two  apparent  improvements 
to  the  imputation  scheme  and  the  results  invite  further  investigation.  A  logical  step 
could  involve  a  more  sophisticated  scheme  that  capitalizes  on  the  hole-boundary 
correlations  between  maximum  S 4  and  <j|4  to  impute  more  realistic  S4  values. 

As  imputation  is  improved,  other  years  and  stations  in  the  dataset  await  perusal. 
Continued  study  and  updating  of  the  information  gathered  by  the  SCINDA  network 
is  encouraged,  along  with  the  generation  and  maintenance  of  statistics  such  as  those 
presented  here.  The  data  should  be  subjected  to  more  robust  statistical  techniques, 
such  as  principal  component  analysis,  to  search  for  potential  relationships  among 
various  stations  and  variables.  The  HD  OP  statistics  in  this  research  were  produced 
using  coarse  almanac  data.  Precise  ephemerides  can  provide  the  means  to  develop 
more  accurate  statistics.  Since  the  HDOP  statistics  generated  were  at  one-minute 
resolution,  it  might  be  instructive  to  calculate  statistics  at  a  15  minute  resolution  for 
each  week,  as  we  did  with  the  availability  data. 
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The  synoptic  climatology  of  outages  could  also  be  examined.  To  our  knowledge, 
no  such  research  has  been  conducted  using  data  from  GPS  receivers.  Rather  than 
attempting  to  fill-in  missing  data,  it  may  be  instructive  to  begin  developing  statistics 
related  to  the  temporal  and  spatial  patterns  of  missingness.  This  should  be  done  not 
only  for  the  types  of  receivers  currently  being  employed  in  the  SCINDA  network,  but 
also  for  typical  receivers  used  by  the  civilian  and  military  community.  Such  research 
may  make  it  easier  to  raise  awareness  of  ionospheric  impacts  on  GPS. 


In  addition  to  work  in  and  among  the  stations  in  the  SCINDA  network,  com¬ 
parisons  can  be  made  between  SCINDA  data  and  remote  sensing  data  collected  by 
platforms  such  as  C/NOFS.  | de  La  Beaujardiere  et  at.  2004  notes  that  C/NOFS 
will  be  monitoring  lightning,  which  has  been  linked  to  “explosive  equatorial  spread 
F”  [  Woodman  and  Kudeki ,  1984  .  Lightning  data  sets  already  exist  and  correlation 
between  lightning  and  equatorial  scintillation  should  be  examined  further. 
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Appendix  A.  Raw  Data  File  Example 


This  appendix  shows  an  example  of  a  raw  data  file  collected  at  Ascension  Island, 
U.K.  during  2001.  These  data  hies  normally  log  an  entire  UTC  day,  but  this 
example  has  been  shortened  and  begins  at  2247  UTC,  or  80280  seconds  and  ends  four 
minutes  later.  PRN  (satellite)  number  6  is  currently  being  scintillated.  The  S4  is 
0.680  at  80820,  rising  as  high  as  0.87  two  minutes  later,  before  being  dropped  for  a 
minute  -  note  PRN  6  is  missing  at  81000,  only  to  return  a  minute  later.  At  81060 
PRN  24  is  now  reporting  an  S4  of  9.00.  Note  in  the  preamble  “9.99=bad  data”. 


This  is  an  error;  “bad  data”  was  actually  assigned  9.00.  See  Section  |3.2|  for  more 
information. 


4;  Number  of  Header  Records 
;  GPS  Computed  Ephemeris  for  1-minute  Data  files 
;  8  Records  per  1-minute  data  record 
-7.94  ;  Station  Latitude 
-14.39  ;  Station  Longitude 
9;  Number  of  Records  per  Measurement  Set  1) 

Yr, dy, sec, nprn, iprn (nprn)  (i2, lx, i3, lx, i6 , lx, i4, etc)  2)  s4 
scintillation  indices  (9.99  =  bad  data)  3)  nprn  repeats  of 
satellite  latitude  (deg)  (f9.3)  4)  nprn  repeats  of  satellite 
altitude  (emr)  (f9.4)  5)  nprn  repeats  of  satellite  longitude  (deg) 
(f9.3)  6)  nprn  repeats  of  azimuth  angle  (deg)  (f9.3)  7)  nprn 
repeats  of  elevation  angle  (deg)  (f9.3)  8)  nprn  repeats  of  300  km 
pp  latitude  (deg)  (f9.3)  9)  nprn  repeats  of  300  km  pp  longitude 
(deg)  (f 9 . 3) 


(A.l) 


1  89  80820  6  4  2  26  7  8  27 

0.680  0.030  0.040  0.030  0.020  0.160 
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10. 

244 

-41. 

176 

-11. 

612 

-45. 

.504 

10. 

228 

30. 

.778 

4. 

192 

4, 

128 

4. 

121 

4. 

.120 

4. 

175 

4. 

.230 

-12. 

214 

8, 

724 

-74. 

462 

-21. 

.075 

2. 

423 

16. 

.826 

6. 

830 

152, 

030 

261. 

172 

187. 

.615 

42. 

932 

36. 

.203 

66. 

152 

40, 

301 

17. 

395 

41. 

.527 

58. 

006 

29. 

.  161 

-6. 

814 

-10, 

543 

-8. 

941 

-10. 

.746 

-6. 

770 

-4. 

.440 

-14. 

254 

-12, 

983 

-21. 

304 

-14. 

.772 

- 

-13. 

295 

-11. 

.827 

89  80880 

6 

4 

2 

26 

7 

8 

27 

0. 

850 

0, 

.030 

0. 

050 

0. 

030 

0. 

010 

0. 

.220 

9. 

837 

-40, 

.862 

-11. 

193 

-45. 

785 

10. 

633 

31. 

.  132 

4. 

191 

4, 

.  128 

4. 

120 

4. 

120 

4. 

175 

4. 

.230 

-12. 

177 

9, 

.008 

-74. 

409 

-20. 

713 

2. 

470 

16. 

.964 

7. 

105 

151, 

.453 

261. 

661 

187. 

128 

42. 

353 

36. 

.000 

66. 

664 

40, 

.438 

17. 

431 

41. 

251 

57. 

594 

28. 

.775 

-6. 

841 

-10, 

.517 

-8. 

882 

-10. 

776 

-6. 

741 

-4. 

.382 

-14. 

252 

-12, 

.963 

-21. 

301 

-14. 

751 

- 

-13. 

290 

-11. 

.804 

89  80940 

6 

4 

2 

26 

7 

8 

27 

0. 

870 

0, 

.020 

0. 

040 

0. 

.010 

0. 

020 

0. 

.200 

9. 

429 

-40, 

.545 

-10. 

775 

-46. 

.064 

11. 

039 

31. 

.485 

4. 

191 

4, 

.  129 

4. 

120 

4. 

.120 

4. 

175 

4. 

.230 

-12. 

141 

9, 

.287 

-74. 

357 

-20. 

.344 

2. 

517 

17. 

.106 

7. 

391 

150, 

.874 

262. 

151 

186. 

.644 

41. 

790 

35. 

.801 

67. 

177 

40, 

.574 

17. 

463 

40. 

.976 

57. 

178 

28. 

.389 

-6. 

869 

-10, 

.491 

-8. 

822 

-10. 

.805 

-6. 

711 

-4. 

.322 

-14. 

250 

-12, 

.944 

-21. 

298 

-14. 

.730 

- 

-13. 

285 

-11. 

.780 

89  81000 

5 

2 

26 

7 

8  27 

0. 

010 

0, 

.040 

0. 

.030 

0. 

.040 

0. 

200 

-40. 

226 

-10, 

.356 

-46. 

.339 

11. 

.444 

31. 

837 

4. 

130 

4, 

.  120 

4. 

.  120 

4. 

.174 

4. 

230 
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9. 

.561 

-74, 

.306 

-19. 

.969 

2. 

.565 

17. 

,251 

150. 

293 

262. 

.641 

186. 

.163 

41. 

.242 

35. 

,606 

40. 

.709 

17, 

.493 

40. 

.701 

56. 

.761 

28. 

,002 

-10. 

465 

-8, 

.763 

-10. 

.834 

-6. 

.681 

-4. 

,262 

-12. 

924 

-21. 

.296 

-14. 

.708 

-13. 

.280  - 

-11. 

.756 

89  81060 

6 

2 

26 

7 

8  27 

24 

0. 

010 

0, 

050 

0. 

.030 

0. 

.040 

0. 

,  170 

9. 

,000 

-39. 

906 

-9. 

937 

-46. 

.611 

11. 

.849 

32. 

,188 

36. 

.601 

4. 

131 

4. 

120 

4. 

.  120 

4. 

.174 

4. 

,230 

4. 

,195 

9. 

829 

-74, 

256 

-19. 

.587 

2. 

.614 

17. 

,399 

-34. 

,870 

149. 

.710 

263, 

133 

185. 

.684 

40. 

.708 

35. 

,416 

337. 

,977 

40. 

845 

17. 

520 

40. 

.427 

56. 

.340 

27. 

,616 

29. 

,518 

-10. 

438 

-8, 

703 

-10. 

.863 

-6. 

.651 

-4. 

,200 

-3. 

,976 

-12. 

905 

-21. 

294 

-14. 

.686 

-13. 

.274  - 

-11. 

.730 

-15. 

,995 
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Appendix  B.  Z-Test  on  Means 

Results  of  Z-Test  on  mean  differences  between  imputed  and  raw  values  for  each  week 
at  a  given  location  for  a  given  year  are  shown  in  Table  [T[  The  first  two  columns  show 
the  percentile  and  S4  level  compared.  The  week  column  refers  to  the  week  in  which 
the  highest  mean  was  recorded  using  the  original  imputation  scheme.  Lower  and 
upper  95%  confidence-interval  limits  for  the  mean  difference  between  imputed  and 
raw  values  are  shown  as  Cl  min  and  Cl  max.  The  Original  Scheme  refers  to  the  re¬ 
sults  obtained  using  the  original  imputation  scheme  while  the  Modified  Scheme  refers 
to  the  results  from  the  imputation  scheme  with  the  addition  of  a  0.3  S4  threshold 
and  the  expansion  of  the  duration  limit  to  20  min.  In  general  the  confidence  inter¬ 
val  limits  increased  for  Anomaly  Crest  stations  (Ascension  and  Parepare)  using  the 
modified  imputation  scheme,  particularly  during  the  solar  maximum.  At  the  stations 
closer  to  the  geomagnetic  equator  (Ancon  and  Marak  Parak),  confidence-interval  lim¬ 
its  decreased  or  remained  the  same.  Although  not  shown,  the  weeks  for  which  the 
highest  mean  was  recorded  using  the  modified  imputation  scheme  tended  to  be  better 
clustered  around  the  equinoxes. 
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Table  1:  Results  of  Z-Test  on  mean  differences  between 
imputed  and  raw  values  for  each  week  at  a  given  location 
for  a  given  year.  The  first  two  columns  show  the  per¬ 
centile  and  S4  level  compared.  The  week  column  refers 
to  the  week  in  which  the  highest  mean  was  recorded  us¬ 
ing  the  original  imputation  scheme.  Lower  and  upper 
95%  confidence-interval  limits  for  the  mean  difference 
between  imputed  and  raw  values  are  shown  as  Cl  min 
and  Cl  max.  The  Original  Scheme  refers  to  the  results 
obtained  using  the  original  imputation  scheme  while  the 
Modified  Scheme  refers  to  the  results  from  the  imputa¬ 
tion  scheme  with  the  addition  of  a  0.3  S4  threshold  and 
the  expansion  of  the  duration  limit  to  20  min. 


Percentile  S4 

Week  Cl  Min  Cl  Max 

Original  Scheme 

Cl  Min  Cl  Max 

Modified  Scheme 

Ascension  2001 

50 

0.3 

47 

0.039 

0.041 

0.051 

0.054 

50 

0.5 

46 

0.051 

0.055 

0.061 

0.065 

50 

0.8 

46 

0.055 

0.06 

0.081 

0.09 

75 

0.3 

45 

0.044 

0.046 

0.051 

0.054 

75 

0.5 

43 

0.065 

0.068 

0.068 

0.071 

75 

0.8 

10 

0.086 

0.093 

0.086 

0.093 

95 

0.3 

10 

0.048 

0.05 

0.048 

0.05 

95 

0.5 

10 

0.05 

0.056 

0.05 

0.056 

95 

0.8 

47 

0.062 

0.064 

0.072 

0.079 

- 1 

Continued  on  next  page 
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Table  1  —  continued  from  previous  page 


Percentile 

S4 

Week 

Cl  Min 

Cl  Max 

Cl  Min 

Cl  Max 

Original  Scheme 

Modified  Scheme 

Ascension  2004 

50 

0.3 

45 

0.06 

0.107 

0.052 

0.082 

50 

0.5 

45 

0.045 

0.066 

0.02 

0.023 

50 

0.8 

3 

0.007 

0.008 

0.026 

0.031 

75 

0.3 

45 

0.099 

0.113 

0.101 

0.126 

75 

0.5 

46 

0.091 

0.109 

0.11 

0.145 

75 

0.8 

46 

0.101 

0.122 

0.101 

0.122 

95 

0.3 

46 

0.081 

0.097 

0.121 

0.279 

95 

0.5 

41 

0.099 

0.148 

0.116 

0.35 

95 

0.8 

45 

0.093 

0.129 

0.103 

0.519 

Ancon  1998 

50 

0.3 

7 

0.002 

0.002 

0.0 

0.0 

50 

0.5 

1 

0.0 

0.0 

0.0 

0.0 

50 

0.8 

1 

0.0 

0.0 

0.0 

0.0 

75 

0.3 

35 

0.022 

0.028 

0.022 

0.028 

75 

0.5 

35 

0.022 

0.028 

0.0 

0.0 

75 

0.8 

35 

0.022 

0.028 

0.0 

0.0 

95 

0.3 

45 

0.051 

0.104 

0.007 

0.007 

95 

0.5 

35 

0.058 

0.11 

0.013 

0.015 

95 

0.8 

35 

0.058 

0.11 

0.013 

0.015 

- 1 

Continued  on  next  page 
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Table  1  —  continued  from  previous  page 


Percentile 

S4 

Week 

Cl  Min 

Cl  Max 

Cl  Min 

Cl  Max 

Original  Scheme 

Modified  Scheme 

Ancon  2001 

50 

0.3 

49 

0.016 

0.017 

0.006 

0.007 

50 

0.5 

51 

0.009 

0.01 

0.005 

0.006 

50 

0.8 

1 

0.0 

0.0 

0.0 

0.0 

75 

0.3 

6 

0.011 

0.012 

0.007 

0.007 

75 

0.5 

46 

0.015 

0.016 

0.017 

0.019 

75 

0.8 

51 

0.01 

0.011 

0.01 

0.011 

95 

0.3 

20 

0.074 

0.124 

0.006 

0.007 

95 

0.5 

20 

0.061 

0.137 

0.011 

0.012 

95 

0.8 

20 

0.061 

0.137 

0.018 

0.02 

Parepare  1998 

50 

0.3 

36 

0.005 

0.006 

0.005 

0.006 

50 

0.5 

40 

0.02 

0.022 

0.02 

0.022 

50 

0.8 

22 

0.006 

0.007 

0.006 

0.007 

75 

0.3 

36 

0.014 

0.015 

0.014 

0.015 

75 

0.5 

36 

0.019 

0.02 

0.019 

0.02 

75 

0.8 

36 

0.022 

0.024 

0.021 

0.023 

95 

0.3 

42 

0.025 

0.027 

0.015 

0.016 

95 

0.5 

46 

0.018 

0.02 

0.011 

0.012 

95 

0.8 

36 

0.034 

0.037 

0.033 

0.035 

- 1 
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Table  1  —  continued  from  previous  page 


Percentile 

S4 

Week 

Cl  Min 

Cl  Max 

Cl  Min 

Cl  Max 

Original  Scheme 

Modified  Scheme 

Parepare  2000 

50 

0.3 

13 

0.029 

0.03 

0.034 

0.036 

50 

0.5 

13 

0.022 

0.023 

0.027 

0.029 

50 

0.8 

13 

0.023 

0.025 

0.035 

0.039 

75 

0.3 

38 

0.024 

0.025 

0.031 

0.033 

75 

0.5 

38 

0.036 

0.038 

0.04 

0.043 

75 

0.8 

13 

0.053 

0.058 

0.068 

0.074 

95 

0.3 

47 

0.024 

0.026 

0.041 

0.056 

95 

0.5 

10 

0.048 

0.052 

0.049 

0.053 

95 

0.8 

13 

0.055 

0.059 

0.068 

0.074 

Marak  Parak  2001 

50 

0.3 

15 

0.005 

0.006 

0.004 

0.004 

50 

0.5 

10 

0.006 

0.007 

0.006 

0.007 

50 

0.8 

1 

0.0 

0.0 

0.0 

0.0 

75 

0.3 

25 

0.005 

0.005 

0.003 

0.003 

75 

0.5 

12 

0.01 

0.011 

0.006 

0.006 

75 

0.8 

39 

0.005 

0.005 

0.005 

0.005 

95 

0.3 

22 

0.011 

0.013 

0.007 

0.007 

95 

0.5 

5 

0.04 

0.043 

0.008 

0.009 

95 

0.8 

10 

0.044 

0.047 

0.032 

0.035 

- 1 
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Table  1  —  continued  from  previous  page 


Percentile  S4 

Week  Cl  Min  Cl  Max 

Original  Scheme 

Cl  Min  Cl  Max 

Modified  Scheme 

Marak  Parak  1998 

50 

0.3 

14 

0.005 

0.005 

0.0 

0.0 

50 

0.5 

16 

0.006 

0.006 

0.0 

0.0 

50 

0.8 

16 

0.006 

0.006 

0.0 

0.0 

75 

0.3 

17 

0.006 

0.006 

0.0 

0.0 

75 

0.5 

17 

0.005 

0.005 

0.0 

0.0 

75 

0.8 

12 

0.007 

0.007 

0.0 

0.0 

95 

0.3 

11 

0.013 

0.013 

0.002 

0.003 

95 
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